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Introduction 


Muddy  coasts  are  found  throughout  the  world  (e.g.  Healy  et  al.,  2002).  Some  of  these  coasts 
are  characterized  by  permanent  extensive  thick  mud  deposits,  covering  part  of  the  hinterland 
(coastal  plain)  as  well,  as  in  the  larger  Mississippi  delta,  including  Louisiana  coastal  system, 
and  the  larger  Amazon  delta,  including  the  Guiana’s.  Elsewhere  the  mud  deposits  are  found 
predominantly  offshore,  as  in  Kerala,  India;  at  some  locations  migrating  to  the  shore 
occasionally,  as  at  Cassino  Beach,  Brazil. 

A  spectacular  phenomenon  occurring  in  muddy  coastal  systems  is  the  damping  of  incoming 
waves,  sometimes  to  a  degree  that  locally  completely  still  water  zones  are  formed  (e.g. 
Gade,  1958;  Kurup,  1972).  For  instance.  Tubman  and  Suhayda  (1976)  carried  out 
measurements  in  1976  in  the  East  Bay  near  the  Mississippi  mouth  at  two  stations  at  20  and 
3.5  m  water  depth  and  3.5  km  apart.  The  bed  consists  of  thick  layers  (several  10  m)  of  soft 
mud  with  a  strength  close  to  the  liquid  limit  (1.75  -  2.36  kPa).  From  a  simple  estimation 
ofthe  wave  damping  to  be  expected  by  bed  friction  only,  they  concluded  that  the  wave 
height  was  dissipated  within  the  soft  mud  by  about  50%  over  this  3.5  km  trajectory. 

Recently,  Sheremet  and  Stone  (2003)  analyzed  five  day  wave  measurements  in  early  2001 
at  two  stations  located  at  the  5  m  isobath  in  the  larger  Atchafalaya  Bay.  One  station  was 
located  above  a  sandy  seabed,  and  the  other  over  a  muddy  seabed.  Though  these  stations 
were  150  km  apart,  the  local  bathymetry,  incoming  wave  conditions  and  local  wind  climate 
were  so  similar  that  a  comparison  of  wave  heights  at  these  stations  is  meaningful.  They 
observed  considerably  lower  wave  heights  (more  than  70%  smaller)  at  the  muddy  site  in 
comparison  to  the  observations  at  the  sandy  site.  This  damping  of  wave  height  is  again 
attributed  to  dissipation  of  wave  energy  in  the  mud  layers.  Moreover,  they  observed  wave 
damping  throughout  the  entire  energy  spectrum,  which  is  unexpected  as  one  would  expect 
that  the  smaller  waves  would  not  be  affected  as  much  by  the  muddy  seabed. 

Wells  and  Kemp  (1986)  measured  wave  heights  in  the  coastal  area  of  Surinam  at  three 
stations  at  7.1,  4.7  and  3.1  m  water  depth  approximately  21.1,  11.7  and  4.3  km  off  the 
coast,  respectively.  This  coastal  area  is  characterized  by  thick  layers  of  Amazon  mud 
(thickness  exceeding  a  few  100  m).  They  observed  88%  and  96%  dissipation  of  the  wave 
energy  traveling  from  the  most  seaward  monitoring  station  towards  the  shore.  It  is  noted 
that  also  in  this  case  dissipation  over  the  entire  energy  spectrum  was  observed,  though  the 
longer  waves  dissipated  most. 

Conceptually,  wave  damping  over  muddy  beds  is  characterized  by  a  sequence  of  events. 
First,  small  elastic  deformations  are  generated  within  the  muddy  seabed  by  the  cyclical 
stresses  induced  by  incoming  waves.  At  a  moment,  these  stresses  may  exceed  the  strength 
of  the  bed,  and  internal  failure  occurs  resulting  in  liquefaction  of  (part  of)  the  seabed.  This 
can  be  a  fairly  rapid  process  (i.e.  of  the  order  of  tens  of  seconds,  up  to  a  few  minutes  at 
most,  e.g.  Foda  and  Zhang  (1994);  Lindenberg  et  al.,  (1989)),  but  this  is  not  always  the 
case.  The  bed  can  then  be  characterized  as  fluid  mud.  The  incoming  waves  now  generate 
internal  waves  at  the  mud-water  interface,  which  are  dissipated  by  internal  friction  within 
the  fluid  mud  layer.  More  sediment  may  liquefy  below  this  fluid  mud  layer,  but  this  is  not 
likely  because  of  two  reasons: 

•  the  fluid  mud  layer  damps  the  waves,  hence  the  stresses  within  the  bed, 

•  earlier  events  (stress  history)  will  limit  the  thickness  of  the  bed  that  is  sensitive  to 
liquefaction. 
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Hence,  one  may  infer  that  the  (maximum)  thickness  of  the  liquefiable  bed  is  a 
characteristic  feature  of  a  specific  coastal  system.  This  thickness  is  determined  by  the  mud 
properties  and  the  wave  climate,  which  together  also  determine  the  important  stress  history. 
Moreover,  erosion  (entrainment)  of  the  mud  layer  may  occur  and  has  to  be  accounted  for 
on  a  somewhat  larger  time  scale. 

The  waves  also  induce  horizontal  stresses  on  the  (fluid)  mud  layer,  known  as  radiation 
stresses.  Hence,  the  fluid  mud  may  move  under  the  influence  of  waves,  provided  that  the 
stresses  exceed  the  strength  of  the  mud.  If  the  mud  cannot  move,  for  instance  when  it  is 
attached  to  the  shore,  some  set-up  may  occur  as  under  normal  wave  conditions;  however,  no 
data  are  available  at  present. 

After  passage  of  waves,  the  liquefied  mud  will  completely  or  partly  regain  its  strength.  As 
the  mud  has  a  very  low  permeability,  liquefaction,  wave  damping,  mud  transport  and 
strength  recovery  take  place  under  more  or  less  undrained  conditions,  and  the  density  of  the 
mud  (water  content)  remains  fairly  constant. 

In  this  report  the  development  of  a  hydro-sedimentological  model  is  proposed,  which 
contains  the  entire  wave-mud  cycle,  i.e.: 

1 .  Liquefaction  of  a  muddy  seabed  by  incoming  waves  during  storm  conditions, 

2.  Damping  of  waves  by  the  liquefied  mud, 

3.  Movement  of  the  liquefied  mud  by  the  waves, 

4.  Strength  recovery  of  the  muddy  seabed  after  passage  of  the  storm. 

Section  2  describes  the  set-up  of  this  model  and  the  rationale  behind  the  various  choices 
that  have  been  made.  The  physical-mathematical  description  of  the  various  modules  and 
their  implementation  are  presented  in  Section  3  and  4.  Section  5  through  7  describe  results 
of  the  model  simulations  on  the  mud  behavior  at  Guyana  coast  and  at  Cassino  Beach.  In 
Section  8  the  results  are  summarized  and  discussed  and  recommendations  for  further 
developments  and  testing  are  given. 
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2  Set  up  the  model 


We  focus  on  the  behavior  of  fairly  old  (i.e.  years,  but  not  centuries  or  more),  non¬ 
consolidating  mud  deposits,  which  liquefy  frequently,  i.e.  of  the  order  of  once  a  year  or 
more  often.  The  strength  of  these  mud  layers  under  quiescent  conditions  is  of  the  order  of 
the  Liquid  Limit,  or  a  bit  larger.  We  do  not  consider  very  old  and  stable  geological 
deposits,  as  encountered  on  the  continental  shelf.  This  implies  that  liquefaction  is  a  fairly 
rapid  process,  i.e.  liquefaction  takes  place  in  several  10  -  100  wave  cycles,  and  not  in  many 
1 ,000  to  1 0,000  characteristic  for  said  geological  deposits. 

We  also  do  not  deal  with  very  recent,  soft  deposits,  so  that  strength  recovery  after  the  wave 
events  takes  place  without  significant  consolidation.  This  implies  that  the  entire  wave-mud 
cycle  takes  place  under  undrained  conditions,  and  that  the  dry  density  (water  content)  of 
the  mud  remains  constant. 

We  further  assume  that  the  mud  deposits  are  more  or  less  homogeneous,  not  containing 
(thick)  layers  of  sand  or  of  other  sediment.  Furthermore,  our  approach  will  be  limited  to 
situations  where  the  thickness  of  the  liquefiable  mud  layers  are  small  in  comparison  to  the 
water  depth,  so  that  the  water  movement  itself  is  not  largely  affected  by  the  behavior  of  the 
mud  layers,  apart  from  a  reduced  interfacial  friction. 

These  assumptions  imply  that  the  wave-mud  cycle  can  be  modeled  with  a  two-layer 
approach,  as  the  interface  between  the  liquefied  mud  layer  and  the  overlying  water  column 
will  be  very  stable  because  of  the  large  density  difference  between  the  two  fluids. 

The  effect  of  viscous  wave  damping  over  fluid  mud  will  be  implemented  in  SWAN.  This 
model  requires  input  of  the  thickness,  extension,  density  and  viscosity  of  the  fluid  mud 
layer,  and  of  course  information  on  bathymetry,  incoming  waves  and  wind  stresses  to 
compute  locally  generated  waves.  The  important  feed-back  between  incoming  waves  and 
the  mud  bed  (liquefaction,  damping,  traction,  consolidation)  has  not  yet  been  accounted  for. 
We  propose  to  model  liquefaction  of  the  fluid  mud,  its  movement  and  the  consolidation 
after  the  wave  events,  and  the  feed-back  with  the  waves  (SWAN)  within  the  DELFT3D-2L 
FLUID  MUD  MODEL.  This  model  is  currently  incorporated  as  a  special  in  the  standard  version 
of  DELFT3D.  A  detailed  description  of  the  model  can  be  found  in  the  User  Manual  (Wang 
and  Winterwerp,  1992)  and  in  Winterwerp  et  al.  (2003).  This  model  consists  of  three  layers 
(e.g.  Section  3.1): 

•  an  upper  layer,  the  water  layer,  which  may  contain  dilute  suspensions  of  mud,  and 
which  is  modeled  in  three  dimensions, 

•  the  depth-averaged  fluid  mud  layer  in  which  the  sediment  concentration  is  assumed  to 
remain  constant, 

•  the  bed  layer  (consisting  of  sand  or  well-consolidated  mud,  which  can  be  eroded  and 
can  store  mud. 
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F/£.  2.1 :  Sediment  balance  in  DELFT3D-2L  FLUID  MUD  MODEL. 

Fig.  2.1  presents  a  sketch  of  the  mass  balance  of  the  2L  model  and  Fig.  2.2  presents  a  sketch 
of  the  momentum  balance,  in  .v-direction  only  for  simplicity. 


A  detailed  description  of  the  existing  DELFT3D-2L  FLUID  MUD  MODEL  is  given  in  Section 
3.1.  Extension  of  the  model  will  be  carried  out  in  small  steps  to  guarantee  testable  sub¬ 
products,  not  only  to  test  the  model  implementations  themselves,  but  also  the 
mathematical-physical  concepts,  as  these  are  new  to  a  large  extent.  The  proposed  phasing 
is  described  in  Table  2.1,  including  reference  to  the  sections  with  further  details.  Column  3 
summarizes  the  material  properties  of  the  liquefied  that  have  to  be  specified  by  the  user:  in 
Phase  1  and  2,  the  thickness  and  extension  #.rvv)  of  the  mud  layer,  its  viscosity  p,  strength 
tb  and  density  p  have  to  be  prescribed.  Thickness  and  extension  of  the  mud  layer  is 
modeled  in  Phase  3  of  the  model  development,  whereas  the  rheological  parameters  of  the 
mud  are  modeled  in  Phase  4. 


Table  2.1 

:  Phasing  of  the  model  set-up. 

Phase 

physical  process 

for  the  given 

parameters 

physical  see 

Sections 

1 

mud-induced  wave  damping 

<W),  p,  tb,  pm 

3.3  &  5 

2 

wave-induced  mud  transport 

S(x,y),  p,  tb,  pm 

3.4  &  6 

3 

wave-induced  liquefaction 

p,  tb,  pm 

3.5  &  7 

4 

strength  recovery  after 
storm 

Pm 

3.6 
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3  Model  development 

3. 1  The  2L  FLUID  MUD  MODEL 

Fig.  3.1  presents  a  sketch  of  the  DELFT3D  2L  FLUID  MUD  MODEL  and  the  definition  of  the 
parameters  used. 


The  three-dimensional  water  movement  and  sediment  transport  in  the  upper  layer, 
containing  dilute  suspensions  of  mud  are  implemented  in  the  2L  FLUID  MUD  model.  We 
assume  a  hydrostatic  pressure  distribution,  simplifying  the  vertical  momentum  equation.  We 
use  Einstein’s  summation  convention,  and  the  continuity  equation  for  water  reads: 


=  0  (/  =  1,2,3) 


(3.1) 


where  Uj(x,y^,t)  is  the  three-dimensional  velocity  vector.  The  vertical  velocity  follows 
from  closure  of  the  three-dimensional  continuity  equation.  The  momentum  equation  reads  (/ 


&/=  U): 


du,  .  du,  f  gh  dp  dZs  d  .  ,  dir 


ph 


■=z, 


ph 


2pdxi 

=  0 


dx,  dxt 


dx. 


(3.2) 


:=Z, 


where  /  =  1,2  (i.e.  x,y  -  see  Fig.  3.1),  e-^y  is  the  permutator,  /  =  Coriolis  parameter,  Z,  is 
either  Zm  if  fluid  mud  is  present,  or  Zh  if  no  fluid  mud  is  present.  Subscript  and  refer 
to  water  surface,  water  -  fluid  mud  interface  and  bed,  respectively.  Further,  v  =  kinematic 
viscosity  and  the  eddy  viscosity  Vj  follows  from  the  standard  turbulence  model  within 
DELFT3D  (for  instance  the  k-E model),  which  is  not  treated  here. 
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The  mass  balance  for  suspended  sediment  reads: 


fic 

—  + 
dt 


fix,  fir. 


(D+rr)f- 

fir. 


=  -D 


"'I  z=Z„ 


(3.3) 


Where  is  the  Kronecker  delta,  D  =  sediment  diffusion  coefficient  and  Tr  is  the  eddy 
diffusivity. 

For  the  mass  balance  for  salt  the  standard  formulation  in  DELFT3D  is  used.  The  equation  of 
state  accounts  for  effects  of  salinity  (and  temperature)  only,  not  for  the  effects  of  suspended 
sediment,  i.e.  p  =  p(S,T). 

The  fluid  mud  is  assumed  to  appear  in  patches,  as  sketched  in  Fig.  1  and  2,  which  can  move 
freely  when  tractive  forces  are  large  enough,  both  along  and  perpendicular  to  the  coast  (N.B. 
this  implies  no  wave-setup  for  the  time  being).  The  behavior  of  the  lower,  fluid  mud  layer 
is  governed  by  the  momentum  equation  and  the  mass  balance  for  sediment  in  the  lower 
layer;  the  latter  simplifies  to  an  equation  for  the  layer  depth,  as  the  sediment  concentration 
in  the  fluid  mud  layer  (bulk  density)  is  assumed  to  be  constant: 


fifi  (  &umJ8 

dt  fix,. 


=  {D„,  ~Em  +  Eh -Wl  )/cm 


(3.4) 


The  momentum  equation  reads: 
fit/ 


dt 


„  of  g Ap 

- eijkfkUn,J+S 


dZ...  dZ  A 


m.j 


'J  fix. 


fix 


fix,  fix. 


d~u.. 


-v 


fix; 


r*-'.+2V^o 


pmSm  pS 

r  m  m  r  m  m 


(3.5) 


The  bed  shear  stress  is  given  by  a  Bingham-kind  rheological  model: 


rA,-  = 


A  2 

TB+—P,nUn 


brelUm,i 


X+bre,^ 


(3.6) 


A  pure  Bingham  model  cannot  be  applied,  as  the  Bingham  stress  would  induce  an  erroneous 
tractive  force  on  the  mud  at  small  flow  velocities.  Therefore,  we  have  added  the  relaxation 
brei,  which  has  a  value  of  about  105  —  106  s/m  to  obtain  finite  values  of  Tz  in  the  rang  of  0 

<  u,„  <  umje\  m/s.  For  instance,  for  brei  =  105  s/m,  Tz>  varies  linearly  from  r7  t//  =0  for 

Km.rel  =  0  to  TZb.eff  =0-99 TZb  at  u„,.rei  =  1  mm/s.  The  stresses  at  the  water-mud  interface 
read: 


r 


Zm.i 


(3.7) 


The  shallow  water  equations  discretisized  in  the  2L  FLUID  MUD  MODEL  are  solved  with  the 
finite  difference  Alternating  Direction  Implicit  (ADI-)method  (Stelling  et  al.,  1986)  on  the 
staggered  grid,  sketched  below.  This  staggering  implies  that  the  layer  thickness  at  the 
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velocity  points  has  to  be  obtained  by  averaging  depth  and  water  level  values  from 
neighboring  points.  This  procedure  was  selected,  as  it  allows  for  internal  super-critical  flow. 
The  thickness  of  the  water  layer  h"  for  the  velocity  point  Uj+aj  for  instance  is  presented  in 
equ.  (3.8)  as  an  example. 


Zs  +  Zs 

Ijt  _  4i j  Vi  j 


zb  +  zh 


(3.8) 


The  flooding  and  drying  procedure  implies  that  the  flow  velocity  is  set  to  zero  when 
Ky<^  where  is  a  user-defined  threshold  value.  Such  an  inactive  velocity  grid  point 

becomes  active  again  when  h'fy  >  2eh . 


A  net  flow  into  a  computational  box  not  necessarily  implies  positive  water  depth  in  the  box 
centre.  Therefore  an  additional  requirement  has  to  be  met:  h?  >  0 ,  where  superscript  C 
refers  to  the  water  level  point.  The  depth  at  a  water  level  point  is  given  by: 


=  max 


(3.9) 


If  this  would  result  in  a  non-positive  water  volume  in  an  active  computational  box,  the 
procedure  at  this  half-time  step  is  repeated,  but  now  with  the  computational  box  set  dry.  The 
flow  velocity  is  set  to  zero  in  the  four  surrounding  velocity  points. 


Zs 


j+!4  o  |  o  |  o 


h 


V 


zm 


zh 

reference  level 


Definition  of  levels  -  see  also  Fig.  3 


)-+-+- 

j-Vl  O  |  O  I  o 

i-Vi  i  i+'/2  i+1 


-  =  //-velocity  component 
|  =  v-velocity  component 
o  =  elevation  of  consolidated  bed  Z/, 

+  =  water  or  fluid  mud  elevation  Zs,  Zm 


The  upwind  scheme  for  super-critical  mud  flow  implies  that  the  mud  thickness  is  defined  as 
a  function  of  the  mud  flow  direction.  If  u l+y  .  >  0  (to  the  right  in  the  sketch  above),  the 

thickness  of  the  fluid  mud  at  this  velocity  point  becomes  =  Zm  -  Zb  ,  whereas  for 

ui+y  j  <  0:  Sui+y  j  =  Zm  -  Zb  .  At  dry  points  the  thickness  of  the  fluid  mud  layer  is 

defined  as  a  function  of  its  slope,  i.e.  tf+XJ  =  Zmu  ~  Zb,  for  Zm,  ,  >  Zm,.1J  and 

d'i+Yuj  =  Z«MJ~ZbM-  f°r  Z'”„  <  Zm  .  Note  that  mass  is  always  conserved  with  this 
procedure. 
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For  the  time  integration  the  ADI-scheme  is  applied.  For  the  one-dimensional  case,  the 
continuity  equation  becomes: 


8>f  -8'j 

At 


+ 


o/  <-/ 

°i+ZJui+HJ  ~  °i-ZJui-ZJ 

Ax 


=  0 


(3.10) 


It  can  be  shown  (Stelling  et  al.,  1998)  that,  for  the  upwind  approach,  positive  definite  water 
levels  are  guaranteed  if  the  Courant  condition  for  velocity  is  met: 


A  tu. 


i+Zj 


Ax 


(3.11) 


3.2  Modeling  water-bed  exchange  processes 


The  exchange  of  sediment  at  the  interface  between  the  consolidated  bed  and  the  dilute 
suspension  is  modeled  as  follows: 

Deposition  from  dilute  suspension: 


Z) 


=  w,c, 


■Z„ 


(3.12) 


Where  Ws  [m/s]  =  settling  velocity  of  the  mud  floes  in  the  dilute  suspension  and  c 
[kg/nr  ]  is  the  suspended  sediment  concentration  close  to  the  bed. 

Erosion  by  dilute  suspension: 

El  I  -r, .)  for  r.>  r 

I  z=Zh  E\l’  c,e  )  b  c.e 


where  r/,  =  bed  shear  stress,  induced  by  currents  and  waves.  The  erosion  parameter  A//. 
[kg/nr/s/Pa]  and  the  critical  shear  stress  for  erosion  rc  e  [Pa]  follow  either  from  user  input, 
or  are  established  through  the  formula  by  Van  Kesteren  (Winterwerp  and  Van  Kesteren, 
2004): 


Mc=k, 


CAoPdn 


io  a 


50 


and  Tce  =0. 163A2(PI)084 


(3.14) 


Where  cv  [m2/s]  =  the  consolidation  coefficient,  (f>s,o  [-]  =  initial  (i.e.  before  swelling) 
volumetric  concentration  of  the  solids  in  the  bed  (solid  fraction),  /}/„  [kg/m3]  =  dry  density 
in  the  bed,  D50  [m]  =  median  size  of  primary  floe  size  in  the  bed,  PI  [%]  =  plasticity  index, 
and  ki  and  kj  are  tuning  coefficients,  if  necessary. 

The  exchange  of  sediment  at  the  interface  between  the  fluid  mud  layer  and  the  consolidated 
bed  and  the  dilute  suspension  is  modeled  as  follows: 
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Deposition  on  fluid  mud: 

D  I  =Wc  7  (3.15) 

m  \z-Zm  s  -=Zm 

Where  Ws  [m/s]  =  settling  velocity  of  the  mud  floes  in  the  dilute  suspension  and  c_=z 

[kg/m3]  is  the  suspended  sediment  concentration  in  the  dilute  suspension  close  to  the  fluid 
mud  layer. 

Consolidation  of  fluid  mud: 

The  consolidation  rate  is  derived  from  Winterwerp  and  Van  Kesteren  (2004),  where  it  is 
shown  that  the  vertical  velocity  ws  of  mud  floes  in  a  consolidating  bed  with  respect  to  a 
Eulerian  reference  frame  is  described  by: 


vv 


—  Ps  Pw 
Pw 


k 

gp»  & 


(3.16) 


Where  <f>  [-]  =  solid  fraction  (=  cjps,  with  c  [  kg/m3  ]  =  mass  concentration),  k  [m/s]  = 
permeability  and  <rl*  [Pa]  =  the  effective  or  skeleton  stress.  We  apply  fractal  theory  to 

establish  k  and  <T„  : 


2  2 

k  =  Kk  ((/>)  i-*f  and  erf*  =  Kp  (V)3-"/  ;  <j)  =  cjps 


(3.17) 


Where  n/  is  the  fractal  dimension,  and  Kk  and  Kp  are  coefficients  to  be  detennined  from 
consolidation  experiments.  Substitution  into  (3.16)  and  averaging  over  the  fluid  mud  thickness 
<!>  yields  (assuming  that  (ft  docs  not,  or  only  slowly  varies  with  time  and  depth): 


r,=^-Jv,d z=k, 


Ps  -  Pv 

pw 


.  KkK„ 

+“ - 


8gPw 


(3.18) 


Typical  values  for  the  various  parameters  amount  to:  rtf=  2.5  -  2.7,  Kk  =  10  n  -  1014  m/s 
and  Kp  =  107  -  1011’  Pa.  This  implies  that  for  c  =  100  -  300  kg/m3,  Vjcm  amounts  to  10‘8  - 
10  111  m/s,  hence  the  rate  of  consolidation  can  be  neglected  in  applying  the  2L  FLUID  MUD 
MODEL  for  the  cases  treated  in  this  report.  This  agrees  with  our  assumptions  that  the  fluid 
mud  behaves  as  an  undrained  fluid. 


Erosion  by  fluid  mud: 

The  same  formula  applied  for  the  dilute  suspension,  i.e.  equ.  (3.13)  can  be  applied  for  the 
erosion  of  the  consolidated  bed  by  the  fluid  mud  layer  as  well.  We  note  that  the  flow 
velocities  of  the  fluid  mud  layer  for  the  cases  considered  in  this  report  are  so  low  in  general 
(often  within  the  laminar  flow  regime),  that  the  erosion  of  the  sub-soil  may  be  negligible  in 
most  cases. 
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Entrainment  of  fluid  mud: 

As  mentioned  above,  the  flow  velocity  of  the  fluid  mud  layer  is  very  small  in  general,  often 
within  the  laminar  flow  regime.  Erosion  of  the  fluid  mud  layer  therefore  takes  place  through 
entrainment  by  the  turbulent  water  movements  over  the  fluid  mud  layer,  and  the  rate  of 
entrainment  reads  (Winterwerp  and  Van  Kesteren,  2004): 


E,„  =u. 


0.5 


5.6  +  Ri , 


(3.19) 


Where  u*  [m/s]  is  the  relevant  shear  velocity  and  Ri*  the  bulk  Richardson  number  defined 
as: 

~  A, 

py* 


The  shear  velocity  u*  has  a  flow-induced  component  and  a  wave-induced  component,  as  a 
turbulent  wave-induced  boundary  layer  may  develop  at  the  interface: 


u:  =u;  f+u;K 

(3.21) 

=  Wl 

(3.22) 

where/  =  friction  coefficient  {f,  ~  4- 10'4),  Ui  =  depth-averaged  flow  velocity  in  upper  layer 
(e.g.  Delft  Hydraulics,  1974),  and  for  the  wave-induced  friction  we  follow  Soulsby  and 
Clarke  (2005): 

Wi,h  /.,=  0.0521  Re-01*’  ,  =  (3 .23 

2  v 

Where  Uorb,i  =  amplitude  wave  orbital  motion  and  A  =  wave  semi-orbital  excursion. 


)gh 


(3.20) 


3.3  Mud-induced  wave-damping 

The  two-layer  fluid  mud  model  by  Gade  (1958)  is  modified  a  bit  and  implemented  in  the 
SWAN-wave  model.  The  SWAN-model  is  a  two-dimensional  fully  spectral  state-of-the-art 
wave  propagation  model  developed  by  Delft  University  of  Technology  (Booij  et  al.,  1999; 
Zijlema  and  Van  der  Westhuijsen,  2005).  Swan  solves  the  balance  equation  for  wave 
action  density  N\c r,  9)  [m2]  per  frequency  cr(=  2iz!T)  and  direction  6,  which  reads: 


D  N' 

D E'  _  S' 

S' 

° h.m 

C' 

.  ^nl 

Dt 

crDt  <7 

<j 

(7 

,  +  _"L 

<J 

(3.24) 
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The  wave  action  density  N'  is  defined  as  N'  =  E'/cr ,  where  E'[a,0)  [nrs]  is  the  wave 
energy  density  and  the  total  wave  energy  is  obtained  by  integration  over  frequency  and 


direction: 

2ji  oo 

E  =  J  Jr d(7d0 

0  0 

(3.25) 

The  root-mean-square  of  the  orbital  velocity  at  the  seabed  uorb  and  the  standard  deviation  of 

the  water  elevation  rj  read: 

u]rb  =  f)  T,  E  d<xd0  and  '  ff  =  E 

o  o  sinh-  {kh}  gp  2  gp 

(3.26) 

and  the  significant  wave  height  Hn  follows  from  E  =  pgH^f 8  [J/m2], 

as  the  first 

moment  and  significant  wave  height  are  almost  identical:  Hi:i  ~  //m0.  Note  that  within 
SWAN  the  energy  is  divided  by  pg  to  yielding  the  energy  density  E  [m2],  which  leads  to  the 
other  units  given  in  these  sections.  In  (3.24),  is  a  source  term  per  wave  frequency 

crand  direction  0  accounting  for  wind  input,  S'h(<J,0)  is  a  sink  term  accounting  for  energy 


dissipation  by  bed  friction,  white  capping  and  depth-induced  wave  braking,  S'hm(cr,0) 
(with  dimension  [m2]  in  SWAN)  is  a  sink  term  accounting  for  viscous  dissipation  in  the  mud 
layer  and  S'n,  (  crj))  accounts  for  non-linear  wave-wave  interactions.  The  first  two  and 

fourth  terms  have  been  implemented  in  the  standard  version  of  SWAN  and  are  not 
elaborated  here.  The  third  term  is  modelled  further  to  the  two-layer  approach  by  Gade 
(1958),  where  we  assume  that  Gade’s  approach  is  applicable  per  wave  direction,  and  that 
the  individual  solutions  per  wave  direction  may  be  added  linearly.  We  note  that  the  2nd 
term  (losses  by  friction)  becomes  much  smaller  when  mud-induced  damping  is  relevant,  as 
the  lower  layer  is  now  fluid  (liquefied  mud). 


The  material  derivative  D/Df  in  equ.  (3.24)  contains  advection  of  energy  in  geographical  jc- 
y  space  (i.e.  the  sum  of  the  group  velocity  Cg  and  (depth)  mean  flow  velocity  U),  and  the 
propagation  velocity  in  spectral  cr-6  space.  If  we  assume  that  U  «  Cg,  only  advection  by 
the  group  velocity  remains. 


Fig.  3.2  shows  the  lay-out  of  the  two-layer  system  and  the  various  definitions:  the  fluid 
mud  layer  depicts  a  bounded  internal  wave.  Free  internal  wave(s)  in  the  fluid  mud  layer  are 
irrelevant  for  the  current  problem,  as  these  are  far  too  long,  and  are  therefore  not  treated 
here.  It  is  important  to  note  that  this  implies  that  the  frequency  a  and  (complex)  wave 
number  k  (=  kr  +  i&,)  are  identical  for  the  upper  water  layer  and  lower  fluid  mud  layer. 
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Fig.  3.2:  Sketch  of  two-layer  fluid  mud  system  and  definitions  (after  Gade ,  1958). 


Gade  (1958)  limits  his  analysis  to  shallow  water  waves  and  he  assumes  that  the  upper  layer 
is  non-viscous,  and  that  the  waves  are  sinusoidal  and  small  compared  to  the  water  depth,  so 
that  vertical  accelerations  can  be  neglected.  Gade’s  model  was  extended  by  Dalrymple  and 
Liu  (1978)  accounting  for  vertical  accelerations  (so  the  model  becomes  applicable  for  short 
waves  as  well)  and  a  viscous  upper  layer;  by  De  Wit  (1995)  accounting  for  vertical 
accelerations;  by  MacPherson  (1980)  accounting  for  vertical  accelerations  and  elastic 
effects  in  the  lower  layer;  and  by  Maa  and  Mehta  (1990)  accounting  for  elastic  effects  in  a 
multi-layered  modelled  fluid  mud  layer.  Here,  we  follow  De  Wit  (1995)  as  his  approach 
matches  the  next  phase  of  our  model  development  best.  The  two-layer  system  is  then 
described  with  the  following  set  of  equations  for  continuity  (equ.  (3.27)  and  (3.29))  and 
conservation  of  momentum  (equ.  (3.28)  and  (3.30))  for  the  upper  (subscript  i)  and  lower 
layer  (subscript  2),  respectively: 


II 

o 

(3.27) 

^l+_L^l  =  0 

dt  px  dx 

(3.28) 

cHv.  1  dp. 

— L+ - —  =  -e 

dt  pt  dx 

du2  dw2  _ 
dx  dz 

(3.29) 

du2  1  dp2  d2u2 

dt  '  pm  dx  ~'m  dz2 

(3.30) 

*1? 

+ 

- 

II 

1 

&Q 

in  which  pi  and  p2  are  the  pressure  in  the  upper  and  lower  layer,  ux  and  z/2  are  the  horizontal 
flow  velocity  in  the  upper  and  lower  layer,  w\  and  w2  are  the  vertical  flow  velocity  in  the 
upper  and  lower  layer  with  thickness  h  and  S,„ ,  respectively,  and  vm  is  the  (constant)  fluid 
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mud  viscosity.  The  horizontal  and  vertical  co-ordinate  are  given  by  x  and  z  (note  that  x  is  in 
the  direction  0  of  the  waves),  and  t  is  time. 

This  set  of  equations  is  solved  with  the  following  boundary  conditions: 

•  at  free  surface  (z  =  h  +  S,„):  p  =  0,  (77,  =  77,  exp  ji(Ax-crr)j)  and  vc,  =  dp Jdt , 

•  at  interface  (z  =  S,„):  dp /dz  =  0 ,  du2/dz  =  0  and  w,  =  w2  =  <3  77,  /dt ,  and 

•  at  consolidated  bed  (z  =  0):  Uj  =  W?  =  0. 

where  <7  is  the  (real)  wave  frequency  and  k  the  complex  wave  number.  Gade  solved  these 
equations  with  the  harmonic  method,  yielding  complex  wave  number  and  amplitudes  of 
layer  thickness  and  velocity  (see  also  De  Wit,  1995).  Gade  presented  results  for  a  uni¬ 
directional,  mono-chromatic  wave  field.  However,  in  SWAN  a  formulation  for  the  energy 
dissipation  S'bm(i 7,9)  is  required.  The  energy  dissipation  was  obtained  by  Gade  by 

integrating  the  work  done  by  the  surface  waves  over  the  wave  period.  As  said,  the  relevant 
wave  height  was  obtained  from  the  harmonic  solution  of  the  set  of  equations  (3.27)  to 
(3.30).  Gade  (1958)  obtained  the  dissipation  per  unit  area  and  unit  time,  averaged  over  a 
wave  cycle  from  the  mean  work  ( fV)  done  by  the  surface  waves  on  the  internal  waves  per 
unit  area: 


Where  p  =  the  pressure  and  n  =  the  normal  vector  to  the  interface,  which  is  correct  for  <£„ 
«  Am.  Extending  Gade’s  analysis1’,  assuming  that  superposition  of  the  solutions  is 
allowed,  the  energy  dissipation  per  wave  frequency  crand  direction  0  reads: 

SL  (£  ■ cr)  =  ocgh0crR^-  sin  {</>'- </>}E'  (3.3 1) 


The  complex  wave  number  ( k  =  kr  +i kj)  follows  from  De  Wit  (1995): 


where  Ag  =  (1  - i)yja/2vm  ,  from  which  kr,  k ,  and  their  modulus  R  +kf  )/< r:) 

and  argument  (phase  angle  between  water  elevation  and  flow  velocity)  </> 
(=2atan  jA:(  /&, .  jj  can  be  computed  |i.e.  ( k/a)'  =/?exp{/^}|.  De  Wit  (1995)  showed 
that  (3.32)  agrees  with  Gade’s  dispersion  relation  for  shallow  water  waves.  For  one-layer 


11  Note  that  Gade’s  paper  contains  a  (typing)error  in  the  elaboration  of  the  energy  losses:  equ.  (3.31) 
should  contain  a  cr. 
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flow  (water  only),  pm  =  p(>  and  Snio  =  0,  and  equ.  (3.32)  reduces  to  the  classical  dispersion 
relation  cr2  =  gk  tanh  kh0 . 

The  relation  between  the  surface  and  internal  waves  follows  from  the  ratio  of  their 
(complex)  amplitudes,  which  can  be  found  from  the  ratio  (  w,  (z  =  Zm  )/vv,  (z  =  Z,  .))  of 

the  vertical  velocity  at  the  water  surface  and  interface,  as  only  bounded  internal  waves  are 
considered.  This  ratio  can  be  obtained  from  De  Wit  (1995)  and  yields: 


*  cr2  1  +  (sinh  {&/?,,}  )’)-g£  sinh  {M0}  cosh  {M0} 

W2  _  __  _  \  v  7  / 

(3.34) 

w,  7,  cr2  cosh{^/70} 

Expanding  (3.34)  around  kho  =  0  with  a  Taylor  series  gives: 

=i-^kh„+Ukh,f  -^K(kK)’ +±-(kh„y +o{(kh,,)’}+ . 

7,  cr  2  6cr  24  '  > 

(3.35) 

which  for  shallow  water  (kho  «  1 )  degenerates  to: 

A  ^  \  2 

— -  =  -r-e'*  =l-^-kh0  =  1  —  gh0  —  =\-ghi)Re"t’  for  kh0  «  1 
w,  rjt  cr  ycrj 

(3.36) 

The  modulus  S  and  argument  (ft  from  equ.  (3.34)  have  been  obtained  with  MAPLE  and 
read2’: 

cr4S2  =  sin  haki  cosh  h0kr  +  gki  cos  h0k:  sinh  h0kr  -  cr2  sin  h()ki  sinh  h0kr ) 

+  [gkr  cos  h()ki  sinh  h0kr  -  gkt  sin  h0ki  cosh  h0kr  -  cr 2  cos  h0kj  cosh  h0kr ) 

+  (3.37 

2  ) 

tan  </>'  =  (gkr  cos  h(lki  sinh  h()kr  -  gkt  sin  h0kj  cosh  h0kr  -  cr2  cos  hnkt  cosh  h0kr ) 
(  gkr  sin  h0k:  cosh  h0kr  +  gki  cos  h0kj  sinh  h0kr  -  cr 2  sin  h0kj  sinh  h()kr ) 

/  (3.38 
) 

Hence,  as  (f)<  \(/>'\  and  ft  <  0  is  always  true,  S'h  m  >  0  (e.g.  Gade,  1958). 

It  is  noted  that  our  analyses  show  some  inconsistencies  in  the  dispersion  relation  and 
the  other  relations,  which  may  be  the  result  of  an  error  in  the  dispersion  relation.  This 
is  being  checked  at  this  moment.  The  dispersion  relation  may  therefore  have  to  be 
altered  in  the  near  future. 


Note  that  Gade's  relations  for  5  and  <j>  read:  R  =  ( k~  +  k  and  <p'  =  2atan  j k  j kr  j . 
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The  mud  properties  Sh,(x,y),  p,„  and  vm  follow  from  measurements  or  the  models  presented 
further  in  this  paper.  Finally,  we  have  added  an  extra  calibration  coefficient  a  in  (3.31)  to 
account  for  possible  non-linear  effects  allowing  a  better  tuning  of  SWAN  to  the 
observations,  if  necessary. 

Gade  (1958)  and  Dalrymple  and  Liu  (1978)  found  that  maximum  wave  attenuation  occurs 
when  the  thickness  of  the  fluid  mud  layer  amounts  to  about  1.3  times  its  boundary  layer 

thickness  Am  ^(r/2vm  j ,  which  agrees  well  with  Gade’s  laboratory  observations  (see 

also  Section  5). 


3.4  Wave-induced  mud  transport 

In  this  section  a  model  for  wave-induced  mud  transport  is  derived.  It  can  be  shown  that  the 
tractive  force  by  wave-induced  interfacial  drag  is  small  compared  to  the  forces  by  the 
radiation  stress.  Only  the  effects  of  radiation  stress  are  therefore  accounted  for. 


Here  we  apply  a  right-handed  co-ordinate 
system  with  positive  z-axis  (x.t-axis) 
positive  upward.  Waves  travel  at  an  angle 
0  with  the  positive  x-axis.  Note  that  X/  is 
not  necessarily  perpendicular  to  the  coast. 


Longuet-Higgins  (1953,  1958)  has  shown  through  a  simple  momentum  balance  that 
momentum  losses  by  bed  friction  in  the  wave  boundary  layer  are  compensated  by  a  tractive 
force,  which  has  the  form  of  the  well-known  Reynolds  stresses.  These  stresses  are 
commonly  referred  to  as  radiation  stresses  and  the  tractive  forces  are  known  as  streaming. 

It  can  be  shown  (see  equ.  (3.41))  that  the  radiation  stress  (i.e.  the  wave  momentum  fluxes)  is 
proportional  to  the  wave  energy  in  the  wave  boundary  layer.  Therefore,  any  change  in  the 
wave  energy,  for  instance  due  to  shoaling,  bed  friction,  viscous  damping,  etc.,  must  induce  a 
change  in  radiation  stress,  hence  tractive  forces.  In  cross-shore  direction  this  accounts  for 
the  well-known  wave  set-up,  and  in  long-shore  direction  in  the  well-known  long-shore 
wave-induced  currents,  which  can  be  quite  large,  in  particular  in  the  surf  zone. 

First,  the  more  general  case  is  treated.  The  hydrostatic  momentum  equation  is  averaged  over 
the  wave  period  (e.g.  Phillips,  1977) 


dua  dua  1  dp 

5.  +UJ  Qx  QaPkfkUa  + 


dr 


ja 


dt 


pdxa  p  dx  / 


_  F  _  1  5SaP 

~ta 


p  dx 


(3.35) 


Thus  for  instance  in  x-direction: 


du  du  du  du 

- ft/ - t- V - +  W - h 

dt  dx  dy  dz 


\_ 

P 


dK^K\ 

dx  dy  j 


(3.36) 
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The  wave-induced  tractive  force  Fa  is  determined  by  the  gradients  in  the  so-called  radiation 
stress  Sap-  Here  a  and  (5  (=  1,2)  correspond  to  the  horizontal  co-ordinates  ,y  and  y  only, 
whereas  j  (=  1,  2,  3)  corresponds  to  the  co-ordinates  x,y,  z.  The  radiation  stress  tensor  Sap  is 
given  by  (e.g.  Mei,  1989;  Dingemans,  1997): 


S  S  1 

XX  AT 

(/7COS2  9  +  n-\l2}E 

(/7Cos#sin#)£ 

s  s 

yx  yy 

(ncos9s\x\9)E 

(rcsin2  9  +  n-\/2^E 

where  n  =  Cg/c  =  kjcr  •  dcr/d/r,  ,  Cg  =  group  velocity,  c  =  wave  celerity.  In  case  of  clear 
water,  we  find  the  common  relation  n  =  /c  =  \/2  + krh/sinh{2krh} ,  from  which,  by 
differentiating  the  normal  XY-component  to  .y,  we  find  for  instance: 


as.  d(  2  .  \  cos  d  6  I 

— —  = — (/jcos  9  +  n-  \  2)E  - - 1  EC  cos 

6x  dx  '  '  c  d x'  g 


d 

krhE  ^ 

cos  9  dPv  d 

1  '  — 

(  k,.hE  ^ 

dx 

v  tanh  {2Ar,./z}  j 

c  dx  dx 

v  tanh  {2 krh)  J 

’  dx 
cos  6 


khE 


tanh  { 2^/7} 


D 


J 


(3.38) 


as  dPa/d. x  +  D  =  0 ,  where  D  =  total  dissipation  rate  per  unit  area.  Hence,  we  find  for  the 
gradient  of  the  radiation  stress  tensor: 


d  |  krhE  I  D cos 9 
dSap _  <5y  jsinh{2&r/t}  J  c 

^xp  D  sin  9 


DcosP 

c 

d  J  krhE  j  Dsm9 
dy  |sinh{2A-,./t}  j  c 


(3.39) 


where  c  =  phase  speed  (celerity),  k,  =  wave  number,  h  =  water  depth,  E  =  wave  energy 
density  =\pgff  fj  =  wave  amplitude,  Hvi  =  significant  wave  height, 

and  9=  angle  between  wave  propagation  and.Y/-axis. 

The  gradients  in  the  normal  component  of  the  radiation  stress  generate  a  wave-induced  set¬ 
up  (set-down)  of  the  mean  water  level,  whereas  the  shear  component  of  the  radiation  stress 
generates  a  tractive  force.  These  terms  should  be  added  to  the  momentum  equation  of  the 
lower  layer  (3.5).  Furthermore,  wave  set-up  due  to  wave  dissipation  in  the  upper  layer 
induces  an  additional  slope  (4th  term  of  (3.5))  -  this  effect  however  is  already  accounted  for 
in  DELFT3D.  To  enable  the  modeling  of  the  effects  of  the  radiation  stresses  on  the 
momentum  balance  within  the  DELFT3D-concept,  we  apply  the  so-called  GLM-approach,  as 
discussed  by  for  instance  Groeneweg  (1999)  and  Walstra  et  al.  (2000).  This  implies  that  the 
momentum  equation  is  re-written  in  the  well-known  Eulerian  way,  but  the  velocities 
computed  refer  to  the  Generalized  Langrangean  Mean  reference  frame,  and  have  to  be 
corrected  for  the  Stokes’  drift  (see  below). 
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Hence,  the  momentum  equation  for  the  lower  layer  then  reads  (compare  with  equ.  (3.5)): 


duL  L  duLm  L  BuLm  ,  dZ .  gAp  d  ,  _  _  ,  d'uL 

dt  m  dx  "  dy  ^  8  dx  pm  dxK  -  J  "  etc2 


(3.40a) 


rm„v 

Tm.Zm.x 

s: 

p  8 

r  m  m 

+ 


™,Zb,x 


p  8 

r  m  m 


=  0 


z=Z, 


dvL  L  dvL  ,  dvLm  ,  dZ  gAp  8 

- —  +  Um  +  Vm  —  -  fui  +  g +  — — 

ot  dx  By  dy  pm  dy 


-7  \ 

(  ~  Zs  )  —  ~ — f" 

\  m  s  )  m  ^^2 


(3.40b) 


rm,y 

Tm.Z„,y 

Pmfim 

p  8 

r  m  m 

m'Zb ,  v 


z=Z.  r  m  m 


p  8 

r  m  m 


=  0 


z=Zh 


where  superscript  L  refers  to  the  GLM  reference  frame.  Note  that  the  force  induced  by  the 
radiation  stresses  BSap/dxp  are  computed  within  SWAN  as  the  FORCE  rm,a  (see  below).  It  is 

further  noted  that  the  wave  field  in  the  real  world  generally  comprises  a  spectrum.  SWAN 
accounts  for  this  spectrum  by  averaging  wave  number  and  wave  period  weighted  for  the 
energy  at  a  specific  frequency.  This  approach  is  applied  here  as  well.  The  mean  real  and 

imaginary  wave  number  kr  and  ki ,  mean  wave  frequency  a ,  mean  wave  direction  9  , 
mean  wave  celerity  c  =  c)  =  <r / kr  and  mean  value  of  the  radiation  stresses  Sap  are  defined 
as  follows: 


/2n  co 


[S'1]  TM0I 
[m1]  WLENMR 
[m1]  WLENMI 


g  —  J  ^cr(9)E'  (<y,9)  dcrdO /  J  j£'(<r,0)  dcrd9 

0  0  /  0  0 
2nx  /2tx  qc 

K-\\K  (cr,Q)E'(<7,d)  dcrdiO  /  derd# 

oo  /  o  o 

In  oo  /2 7i  x 

(<t,0)E'(<j,0)  dcrdO  /  \\e'(cj,0)  dcrdO 

0  0  /  0  0 

)2tix  /2tix 

J  j£,'(<r,^)sin6>dcrd^  /  |  ^E' (cr,0)cosO  dcrdO  l  [rad]  dir 

oo  /  o  o 

271  X 

S„  =  pg  |  j*(/?cos2  0  +  n-\/2^E'(cr,0)  dodO 

0  0 

2n  x 

SXy  =  Syx  =  pg  J/(  ns'm0cos0)E'(a,0)  d<rd<9 

0  0 

2n  x 

Svv  =  pg  J  J(wsin:  0  +  n-\/2^E'(cr,0)  derd# 


(3.4*) 


r  =  — 


0  0 

as.  as,, 


[j/m:  J  computed  in  swan 


dx  dy 


and  r,  =  — 


as  as 

at  n 

dx  By 


[Pa]  FORCE 
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Where  n(cr,0)  =  kr  (<r,#)/cr(#)-dcr(#)/cl£,  (cr,0).  See  Section  3.3  for  definitions. 

2jt  x 

Note  that  the  radiation  stress  Sxx  =  pmg  J  |(a7COS2  6  +  n- 1/2)£,'|  (cr,#)  dcrdO ,  etc.  in 

0  0 

the  lower  layer  is  not  computed  in  SWAN  and  should  be  obtained  from  the  ratio  between  the 
total  wave  energy  in  the  fluid  mud  layer  Em  (c.q.  E'm )  and  the  total  energy  of  the  surface 

waves  E  (c.q.  E ') ,  using  relation  (3.34).  The  same  holds  for  the  tractive  force: 


E,„=^¥E  ,  hence  rma=^?=-ra  (force) 

/w  p,n , 


(3.42) 


The  actual  Eulerian  velocity  of  the  fluid  mud  layer  follows  from: 


Mm.,  =  UL  =  *4/  -  ul,  =  “Lj  K,x  =  K  cos  6 ,  k  =  kr  sin  9 


Pj  <? 


(3.43) 


E  S 

where  superscript  and  refer  to  Eulerian  reference  frame  and  Stokes’  drift,  respectively. 
The  bed  shear  stress  is  related  with  the  Bingham  strength  Tb  through  a  rheological  model 
(e.g.  Section  3.6): 


m.Zh,i 


A,  /  l  \“ 

^+-A„  K,) 


KelUi,i 


1  +  b. 


re l 


(3.44) 


and  the  stresses  at  the  water-mud  interface,  assuming  ui  »  ufm  : 


A  •>  u . 
o  -in 


(3.45) 


3.5  A  criterion  for  liquefaction  of  the  mud  bed 

In  this  section  we  propose  a  method  to  determine  whether  a  mud  bed  liquefies  and,  if  so,  to 
compute  the  thickness  and  the  extension  of  the  liquefied  mud  layer.  This  method  should 
include  the  liquefaction  of  the  muddy  seabed  by  waves  and  the  strength  recovering  after  the 
waves  in  calm  weather  conditions.  Moreover,  the  effects  of  stress  history  should  be 
accounted  for.  A  model  for  strength  recovery  is  proposed  in  Section  3.6. 

At  present,  no  complete  theory  exists.  There  is  some  theoretical  work,  based  on  the  Biot 
model  (Sumer  and  Fredsoe,  2002;  McDougal,  et  al.,  1989)  and  some  experimental  work 
(Alba  et  al.,  1976)  published,  but  none  of  the  studies  are  complete  and  applicable  to  the  2L 
modeling  proposed  here.  Therefore,  we  propose  a  new  method,  which  is  partly  based  on 
earlier  work  carried  out  at  Delft  Hydraulics  (Van  Kesteren  and  Comelisse,  1995). 

It  is  illustrative  to  discuss  the  stress/strength  evolution  in  the  bed  first  in  a  qualitative  way. 
Consider  a  flat  seabed  of  well-consolidated  mud.  The  vertical  skeleton  (effective)  stress  at  a 
certain  depth  is  determined  by  the  density  and  packing  of  the  sediment  and  that  depth.  The 
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lateral  skeleton  stress  follows  from  the  coefficient  of  lateral  stress  (or  coefficient  of  earth 
pressure  at  rest)  Ko: 

c-12  =  (3.42) 


where  07  ( =  J/T^.gdz ,  with  pj,y(z)  =  dry  density  of  sediment  and  z  =  depth  into  the  soil)  is 

the  vertical  principal  skeleton  stress  and  07  and  0>  are  the  two  horizontal  principal  stress 
components.  The  coefficient  of  lateral  stress  Ko  is  related  to  the  angle  of  internal  friction  <p 
through: 


V 


\-v 


l-sin(^) 


(3.43) 


where  v  is  the  Poison  ratio.  The  stress  conditions  in  the  soil  can  be  visualized  in  a  cr-r 
diagram,  the  so-called  Mohr  diagram,  as  depicted  in  Fig.  3.3,  showing  two  circles  of  Mohr 
at  two  different  depths  in  the  seabed,  where  a  =  normal  stress  and  r  is  shear  stress.  For 
simplicity,  we  draw  the  circles  in  one  plane  only.  Note  that  we  depict  the  skeleton 
(effective)  stresses.  For  water,  Ko  =  1  and  07  (=  07)  =  07  and  the  Mohr  circles  degenerate  to 
one  point  on  the  cr-axis;  in  case  <7  represents  skeleton  stresses  only,  G\  (=  07)  =  Oj  =  0,  of 
course.  In  case  of  very  fresh  mud  deposits,  =  0,  and  <T\  (=  07)  =  0. 


Next,  let  us  depict  the  stresses  in  the  seabed  induced  by  regular  waves  with  wave  height  Hs 
and  wave  number  k  and  wave  frequency  cy  in  „v-direction.  Then,  the  wave-induced  water 
pressure  (we  can  omit  the  effects  of  the  wave-induced  bed-shear  stresses,  which  hardly 
generate  stresses  within  the  bed)  on  the  seabed  reads  (Flamant  -  see  Yamamoto  et  al.,  1978; 
Sumer  and  Fredsoe,  2002): 

pgH  cos  (cot  -kx-i//) 

P"  ~  2  cosh  {it/?}  (3-44 

Where  h  the  local  water  depth,  and  cp  is  phase  angle.  The  amplitude  of  the  stresses  pw  in  an 
elastic  bed  reads: 
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pgH  F  Glr 

P"  2  cosh  [kh]  F-  tanh  { A-/?}  pwar  (1  - 1») 

where  G  is  the  elastic  shear  modulus.  For  the  time  being  we  assume  non-elastic  bed 
properties  (G  =  F  =  oo),  e.g.  Sumer  and  Fredsoe  (2002),  so  that  the  stress  amplitude 
becomes  pw  =  pgH/2 cosh  {kh},  and  the  total  skeleton  stresses  within  the  bed  read 
(where  superscript w  refers  to  wave-induced  stresses): 

°zz  =  <  +  \Pdngfc'  =  A,  ( 1  +  kz )  e'fc  cos  ( cot  —  yf )  +  \pJng<h' 

Vxx  =  <  +  Ko  \pjrygte'  =  pw  {\-kz)e-kz  cos((ot-y/)  +  K0  jp^gdz' 

<rt.  =  =  pwkze~b  sin  (cot  -  ip ) 


We  have  also  added  the  effect  of  the  isotropic  stresses  at  depth  z  by  the  own  weight  of  the 
sediment,  with  Pdn(z)  =  dry  density  of  the  sediment  bed.  Perpendicular  to  the  waves  in  the 
horizontal  plane  (v-direction),  the  normal  stress  is  found  from: 


<7. ...  = 


_  “  ,  _ »» 
<T  +  <T__ 


+  K0  jpjngfc’  =  PHe~k:  cos  {(Ot  -  ip) +K0\  pJngdz' 


(3.47) 


If  the  wave-induced  stresses  are  large  enough,  bonds  between  the  floes  in  the  bed  will  be 
broken  and  the  effective  stress  decreases  and  a  water  overpressure  is  building  up.  This  is 
sketched  in  Fig.  3.4,  showing  that  the  original  Mohr  circle  of  Fig.  3.3  (dotted  line)  shifts  to 
the  left,  as  Oj  and  Ko  become  smaller  (hence  ay  decreases  more  than  cr?).  We  note  that  the 
circle  becomes  larger  as  the  shear  stress  increases,  which  is  due  to  the  fact  that  less  bonds 
have  to  carry  the  same  stress.  When  the  mud  lies  on  a  slope,  partial  liquefaction  may  be 
sufficient  to  destabilize  the  bed  (mud  flow). 

Because  of  the  breakage  of  bonds,  also  the  cohesion  and  angle  of  repose  (p  (and  thus  Ko ,  e.g. 
equ.  (3.43))  decrease,  resulting  in  a  translation  and  rotation  of  the  Mohr-Coulomb  failure 
envelope,  as  sketched  in  Fig.  3.4  by  the  dashed  line.  Note  that  we  assume  undrained 
conditions,  which  is  generally  true  for  the  conditions  at  which  the  2L  fluid  mud  model  is 
applicable  (e.g.  Winterwerp  and  Van  Kesteren,  2004). 

The  wave-induced  stresses,  given  in  equ.  (3.44)  should  be  added  as  well  to  the  stress 
diagram  of  Fig.  3.4.  This  is  done  for  one  arbitrary  time,  indicated  by  the  solid  circle  in  Fig. 
3.4  -  basically  a  family  of  circles  is  found,  depending  on  the  phase  of  the  wave.  When  a 
critical  number  of  bonds  is  broken,  the  mud  bed  fails  and  becomes  fluid,  and  the  circles  of 
Mohr  become  smaller  and  smaller  as  the  fluid  mud  is  no  longer  able  to  sustain  considerable 
shear  stresses  (we  get  a  non-Newtonian  fluid);  the  circles  may  even  degenerate  to  a  line  on 
the  cr-axis  around  the  origin  of  the  cr-r  diagram  if  the  fluid  cannot  sustain  any  shear  stresses. 
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After  passage  of  the  waves,  the  mud  can  regain  strength  through  a  kind  of  consolidation 
during  which  the  bonds  between  the  floes  are  restored,  under  quite  conditions.  This  is  depicted 
schematically  in  Fig.  3.5,  showing  that  the  stresses,  angle  of  repose  and  cohesion  are  restored 
(e.g.  Section  3.6). 

As  said,  the  liquefaction  process  was  considered  to  take  place  under  undrained  conditions,  so 
no  water  was  expelled.  However,  in  practice  this  may  not  be  the  case  entirely,  as  a  result  of 
which  the  sediment’s  bulk  density  may  increase  a  bit.  This  implies  that  when  all  bonds  are 
restored  again,  the  bed  may  become  stronger  than  before  the  liquefaction  process.  This  is 
quantified  by  the  so-called  Over  Consolidation  Ration  (OCR),  which  can  be  measured  from 
oedometer  tests  (Lambe  and  Whitman,  1 979). 


In  theory,  the  description  above  can  be  used  to  quantify  the  liquefaction  behavior  or 
liquefaction  potential  of  a  seabed  if  we  know  how  the  skeleton  stresses  decrease  with 
number  of  wave  cycles.  The  latter  can  be  achieved  through  cyclical  roto-viscometer  tests  or 
cyclical  triaxial  tests.  However,  it  is  clear  that  the  use  of  Mohr-circles  is  not  very  practical, 
as  they  have  to  be  drawn  for  each  depth  in  the  seabed  and  for  each  phase  in  the  wave  cycle. 
Therefore,  it  is  more  common  to  visualize  the  stresses  in  the  so-called  p-q  diagram,  where  p 
is  the  isotropic  stress  and  q  is  the  deviatoric  stress.  A  cyclical  load  then  yields  the  stress  path 
depicted  in  Fig.  3.6,  where  the  stresses  follow  typical  butterfly  paths.  When  the  stresses 
exceed  the  failure  envelope,  liquefaction  can  be  expected. 
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Fig.  3.6:  Variation  of  shear  stress  and  water  pressure  during  cyclical  shear  tests  on  Haringvliet  mud  at  various 
shear  angle  amplitudes  (Van  Kesteren  and  Cornelisse,  1995). 

In  our  model  we  make  use  of  the  observation  that  with  the  breakage  of  bonds,  the  water 
overpressure  rises  with  number  of  wave  cycles.  Full  liquefaction  is  found  after  N(  cycles, 
where  N(  can  vary  between  a  few  10  cycles  to  many  thousands  (and  more),  depending  on 
the  wave  characteristics  and  the  strength  (history)  of  the  bed,  e.g.  Sumer  and  Fredsoe 
(2002).  The  results  of  experiments  on  granular  material  (sand)  by  Alba  et  al.  (1976)  are 
depicted  in  Fig.  3.7  as  an  example,  showing  the  so-called  Cyclical  Stress  Ratio  (CSR) 
<7(z)//?o(z),  i.e.  the  ratio  of  the  amplitude  of  the  stress  deviator  in  the  seabed  and  the 

initial  skeleton  (effective)  stress  prior  to  liquefaction  as  a  function  of  the  number  of  wave 
cycles  and  initial  packing. 
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Van  Kesteren  (2000)  carried  out  cyclical  triaxial  tests  on  sediments  from  the  Orange  River,  the 
results  of  which  are  shown  in  Fig.  3.8.  The  results  suggest  a  power  law  relation  between  the 
development  of  the  deviatoric  stresses  and  the  number  of  wave  cycles.  Note  that  liquefaction 
of  the  well-consolidated  clayey-silt  sample  is  only  observed  after  about  60,000  cycles.  As  the 
cycle  period  was  set  at  15  s,  this  would  imply  liquefaction  after  10  days  only.  In  contrast,  the 
sandy  sample  liquefies  after  about  10  cycles,  which  is  in  agreement  with  the  observations  of 
Fig.  3.8. 


These  results  suggest  that  we  can  establish  the  following  criterion  for  failure  of  the  bed,  i.e. 
liquefaction  occurs  when: 


1 


f  N  Y 

q>qcr=Asap0  -j-  or 

J  Asa  p0 


rNV 


>i 


(3.48) 


Where  N  is  the  number  of  actual  wave  cycles,  Nf  is  the  number  of  cycles  at  which  liquefaction 
occurs,  and  a  and  /? have  to  be  determined  empirically  from  cyclical  triaxial  tests  at  a  number 
of  stress  conditions.  The  clayey-silt  sample  of  Fig.  10  would  yield:  a  =  0.6  and  /?  =  -0.24  (e.g. 
Van  Kesteren,  2000).  The  structure  parameter  As  is  discussed  in  Section  3.6;  As  varies  between 
0  and  1 . 


The  isotropic  and  deviatoric  stress  follow  from  (e.g.  Winterwerp  and  van  Kesteren,  2004): 


P  =  /i/3 

^  =  V7.2-3/2 


(3.49) 
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In  which  the  stress  invariants  /,  are  defined  as: 

A  =<Txv+Cr.n+^ 

I 7  =  <T„<TW  +<rn.<7„  +  <7__<7rr  -<Tv.(T_v 


(3.50) 


xv  xz  zx 


as  crv,  =  cr,-  =  0.  After  some  algebra: 
a  =  .  (j2  +  cr  +  <r2  —  cr  a  —cr  a  -cr  cr  +3 cr2 

1  xx  yv  : :  xv  vy  \y  zr  zr  xv  .xr 


(3.51) 


And  upon  substitution  of  the  wave-induced  stresses  into  (3.46)  we  find  for  the  total 
deviatoric  stress  q  and  its  amplitude  q ,  including  the  isotropic  part: 


=  V3 


+— 

3 


p2wk2z2e-2lc  +  pKkzeh  cos  (rtf  -<p)(l-  \pdngdz 

,  V  A'/2 

(l-^o  )\P*8*' 


(3.51a) 


a  =  73 


+  — 
3 


plk2z2e~2k:  +  pHkze~h  (1  -  K0  )|^  fPdngdz' 
\2  V7~ 


+ 


(\-K»)JpJr,gdz 


(3.51b) 


Note  that  the  dry  bed  density  is  kept  constant,  so  that  pj,y  (x,y^,t)  =  p,„  =  const.  The  initial 
effective  (isotropic)  stress  po  prior  to  liquefaction  follows  from  substitution  of  (3.41)  into 
(3.49): 


Po  =  +  °yy  +  <r=  )0  =  2A”  +  1  \p*ygfc' 


(3.53) 


Hence,  we  have  derived  simple  explicit  formula  for  the  deviatoric  and  isotropic  stress 
varying  with  depth  z  into  the  seabed  as  a  function  of  wave  characteristics;  at  depth  z  =  zcn 
the  inequality  of  relation  (3.48)  is  met,  and  the  bed  fails.  We  assume  that  all  sediment  above 
Zaii  becomes  liquefied  and  behaves  as  a  fluid.  Hence,  zan  defines  the  thickness  S  of  the 
fluid  mud  layer  in  the  2L  model. 

The  stress  history  of  the  bed,  which  can  be  quantified  with  the  Over  Consolidation  Ratio 
(OCR)  can  be  modeled  through  a  (equ.  (3.48)).  However,  in  the  cases  we  are  interested,  at 
which  the  fluid  mud  layer  liquefies  and  regains  strength  frequently,  the  OCR  will  be  close  to 
unity  in  general. 

In  case  of  irregular  waves,  a  preliminary  inspection  of  the  wave  spectrum  is  necessary  to 
assess  which  part  of  the  spectrum  is  relevant  for  the  liquefaction  of  the  mud  bed.  This  is  in 
particular  the  case  when  the  spectrum  contains  two  peaks,  with  a  first  peak  generated  by 
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swell  and  a  second  peak  generated  by  local  wind.  In  general,  the  longer  waves  will  be  the 
more  relevant  for  the  liquefaction  of  the  waves.  It  is  not  obvious  which  wave  conditions  are 
critical  for  liquefaction  -  this  may  vary  from  site  to  site,  and  possibly  even  at  a  particular 
site.  The  larger  waves  generate  the  largest  stresses,  but  do  not  occur  frequently.  Smaller 
waves  generate  smaller  stresses,  but  occur  more  frequently.  Therefore,  we  anticipate  that  the 
mean  wave  height  and  corresponding  wavelength  and  period  are  the  dominant  waves  with 
respect  to  liquefaction. 

As  an  example,  the  CSR  is  plotted  in  Fig.  3.9  for  various  wave  characteristics.  We  have  also 
plotted  a  hypothetical  critical  stress  ratio  (q/p)  ,  assuming  that  all  the  mud  above  the 

level  where  the  actual  stress  level  exceeds  said  critical  value  becomes  liquefied.  As  we 
assume  that  liquefaction  takes  place  in  a  time  period  which  is  short  compared  to  the  overall 
storm  period,  we  are  not  concerned  with  the  time  necessary  for  liquefaction,  i.e.  we  do  not 
take  into  account  the  number  of  wave  cycles.  The  critical  stress  ratio  (q/p)  has  to  be 
assessed  experimentally,  i.e.  through  cyclical  triaxial  tests. 


q/p  H 


Fig.  3.9:  Ratio  of  deviatoric  stress  and  initial  effective  stress  (Cyclic  Stress  Ratio)  as  a  function  of  wave  parameters 

and  water  depth. 


3.6  Strength  evolution  and  the  rheological  model 

The  non-Newtonian  behavior  of  the  mud  is  described  with  a  Bingham  model  in  the  2L  FLUID 
MUD  MODEL  (e.g.  equ.  (3.5)  and  (3.6)).  Therefore  we  propose  to  extend  the  Bingham  model 
to  account  for  the  effects  of  waves  and  thixotropy.  The  general  formulation  of  the  Bingham 
model  is  given  by  Malvern  (1969)  and  reads  for  stresses/strains  in  the  vertical  plane: 


u  = 


2//  + 


>7 


D,  ,  i  =  x,y  ;  T/:  >  rB 


(3.54) 
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where  T'r  =  deviatoric  stress  tensor,  //  =  dynamic  viscosity,  xB  =  Bingham  yield  strength,  D,: 
=  rate  of  deformation  tensor,  and  II/j  =  second  invariant  of  Di:;  we  have  assumed 
T'v  « Dxy  » 0 .  Next,  we  assume  that  Di:  consists  of  a  contribution  by  flow-induced 

f 

deformation  (superscript  )  and  wave-induced  deformation  (superscript  "),  and  that  these 
effects  can  be  superposed  linearly:  Dr  =  D  '_  +  D" . 

Furthermore,  we  introduce  a  structure  parameter  to  account  for  thixotropic  effects,  i.e.  the 
muddy  seabed  experiences  repetitive  liquefaction  under  wave  action  and  strength  recovery 
under  calm  weather  conditions.  Our  approach  is  based  on  Moore  ( 1 959),  who  proposed  a 
simple  model  to  describe  such  thixotropic  effects  through  a  structure  parameter  As  (see  also 
Coussot,  1997).  This  model  can  be  extended  (Winterwerp,  1993,  see  also  Toorman,  1997) 
by  incorporating  this  structure  parameter  also  in  the  yield  strength,  which  reads  for  stresses 
in  one  vertical  plane  only: 

r  =  zVoo  +  Aj'b  +(jux  +  Asju')y  (3.55) 

where  x'B  =  yield  (Bighmam)  strength  parameter,  xB,x  =  residual  yield  strength  for  fully 

deflocculated  sediment,//*  =  viscosity  for  fully  deflocculated  sediment,  //'=  a  viscosity 
parameter,  and  y  =  shear  rate.  The  structure  parameter  As  =  1  for  fully  flocculated 

sediments  (all  bonds  active)  and  As  =  0  for  fully  deflocculated  sediment  (all  bonds 
destroyed).  This  behavior  is  depicted  schematically  in  Fig.  3.10. 


shear  stress  r  or  T-.  [Pa] 


shear  rate  y  or  Dy  [s1] 


Fig.  3. 10:  Bingham  model  exle,ided  with  general  Mohr-concepl  account ing  f or  thixotropy 


Combining  these  ingredients  yields: 

-/)-!  (3.56) 


in  which  the  effective  viscosity  /ueff  is  given  by: 


/v  = 


rB.;+Asr'B 

n/iC 


+  /4,  +  As/j' 


(3.57) 
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The  flow-  and  wave-induced  rates  of  deformation  are  then  given  by: 


611 


D!  and  £T=-^- 


(3.60) 


where  q:  is  the  wave-induced  deviatoric  stress  component  in  x-  or  ^-direction  (e.g.  equ. 
(3.46),  (3.47)  and  (3.51b)). 

By  similar  arguments,  we  scale  the  soil  mechanical  parameters  with  the  structure  parameter 
as  well,  so  that  as  well  the  cohesion  c  and  the  angle  of  internal  friction  ^become: 


c  =  Asc0  and  <f>  =  As(p()  hence  K0  =  l-sin(A^>0)  (3.59) 

This  thixotropic  model  can  be  used  to  establish  the  viscosity  and  stress  terms  in  the 
momentum  equation  of  the  2L  FLUID  MUD  MODEL,  e.g.  equ.  (3.5).  Generally,  the  fluid  mud 
flow  will  be  laminar.  This  implies  that  the  friction  coefficient  A  (equ.  (3.6)  and  (3.7))  scales 
linearly  with  viscosity,  i.e.  we  assume: 

A  =  k  /4+V''  (360) 

<*>/("«  )2+(v») 

In  which  kx  is  a  coefficient  to  be  tuned  (from  laminar  flow  theory,  k\  =  64,  so  we  anticipate 
10  <  kx  <  100).  Furthermore,  the  viscosity  and  Bingham  strength  in  equ.  (3.5)  and  (3.6) 
read: 


v  =  (m~  +Kv')!pm  and  rB  =  rBx  +Ast'b  (3.61) 

Note  that  the  viscosity  v  is  also  used  in  the  wave-damping  module  in  SWAN  -  see  Section 
3.3. 

The  rate  of  flocculation  is  given  by  Moore  (1959): 

^  =  0(l-A,)-M,7n7  (3.6!) 

where  a  and  b  are  flocculation  and  floe  break-up  parameters,  to  be  determined  empirically. 
Toorman  (1997)  found  from  roto-viscometer  tests  ^11^,  =  values  of  a  =  Oj  10 2}  s'1  and 

b  =  O { 1  O'4} .  This  implies  a  very  rapid  response  of  the  seabed  to  changes  in  forcing 
conditions.  This  will  have  to  be  verified  for  the  conditions  under  consideration,  but  most 
likely  the  time  scales  will  be  considerably  larger. 


WL  |  Delft  Hydraulics 


3-23 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud- induced  wave  transport  &  wave- induced  mud 
transport 


Implementation  in  DELFT3D 


Finally,  we  have  to  elaborate  (3.54)  for  the  2L  FLUID  MUD  MODEL.  This  implies  that  we  have 
to  integrate  the  wave-induced  deviatoric  stress  q  over  the  mud  layer  thickness  <$ 

Q  =  Jt/dr .  Then  IId,  both  in  equ.  (3.57)  and  (3.61)  becomes: 

s 


nD  = 


fu.  V 


+ 


/v 


(3.63) 


It  is  noted  that  the  denominator  of  (3.57)  becomes  zero  when  Q  =  it  =  v  ,  so  this  will  have 
to  be  implemented  accurately  in  the  numerical  code. 

The  material  parameters  a  and  b  have  to  established  from  rheological  experiments,  for 
instance  through  cyclical  triaxial  tests.  We  note  that  from  small  scale  experiments,  the  time- 
scale  for  recovery  appears  to  amount  of  the  order  of  (a  fraction)  of  seconds.  This  seems 
unrealistically  fast.  However,  from  experience  with  larger  samples,  strength  recovery  is 
expected  to  take  place  within  a  period  of  the  order  of  hours  up  to  a  day  at  most. 
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4  Implementation  into  the  2L  fluid  mud  model 


The  cycle  of  liquefaction,  traction  and  strength  recovery  is  sketched  in  Fig.  4.1a  -  d.  At  /  = 
0  (i.e.  initial  model  conditions)  a  maximum  liquefiable  mud  layer  is  defined  by  the  user 
with  spatial  varying  thickness  ^nax(x,>';/=0).  This  mud  layer  has  the  properties  of  the  fully 
flocculated  material,  i.e.  As  =  1  (e.g.  Fig.  4.1a).  Then,  when  waves  arrive  at  the  relevant 
domain,  all  or  part  of  the  liquefiable  mud  gets  liquefied,  and  the  liquefied  mud  with 
extension  ^T,(.v,jy;/=/|)  attains  the  properties  of  the  fully  deflocculated  material,  i.e.  As  =  0 
(e.g.  Fig.  4.1b.).  Next,  the  liquefied  mud  is  transported  by  the  wave  stresses,  as  sketched  in 
Fig.  4.  lc  (/  =  h).  Then,  strength  is  recovered.  As  =  1  again,  and  a  new  spatial  distribution  of 
liquefiable  mud  is  obtained  (e.g.  Fig.  4. Id,  t^ty). 


Fig.  4.2  presents  the  flow  scheme  of  the  new  2L  FLUID  MUD  MODEL. 


WL  |  Delft  Hydraulics 


4-  1 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud- induced  wave  transport  &  wave- induced  mud 
transport 


Implementation  in  DELFT3D 


Fig.  4.2:  Flow  scheme  o£2L  FLUID  MUD  MODEL 


Finally,  we  repeat  Table  2.1  describing  the  phasing  of  the  model  development. 


Table  4.1 

:  Phasing  of  the  model  set-up. 

Phase 

physical  process 

for  the  given 

parameters 

physical  see 

Sections 

1 

mud-induced  wave  damping 

%X,y),  p,  TB,  pm 

3.3  &  5 

2 

wave- induced  mud  transport 

%X,y),  P,  TB ,  Pm 

3.4  &  6 

3 

wave-induced  liquefaction 

Pi  FBi  pm 

3.5  &  7 

4 

strength  recovery  after 
storm 

Pm 

3.6 
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5  Application  mud-induced  wave  damping 


5. 1  Validation  against  laboratory  experiments 

De  Wit  (1995)  carried  out  wave  damping  experiments  in  a  40  m  long  wave  and  flow  flume, 
with  a  width  and  depth  of  0.8  m.  A  false  floor  was  mounted  on  the  flume  bottom  to  create  a 
0.2  m  deep  test  section  of  8  m  length  in  which  a  dense  mud  suspension  could  be  placed. 
Regular  waves  were  generated  with  a  wave  paddle.  During  the  experiments  wave  height 
and  period  were  measured,  together  with  vertical  profiles  of  the  flow  velocity  and 
suspended  sediment  concentration  at  a  number  of  places,  and  the  total  and  pore  water 
pressure  at  four  heights  within  the  dense  mud  suspension. 

Dense  mud  beds  were  formed  in  the  test  section  through  consolidation  of  more  dilute  mud 
slurries.  In  this  report  we  use  De  Wit’s  experiment  III,  test  3,  carried  out  with  China  clay 
(major  mineral  component  consists  of  kaolinite)  to  validate  the  wave  damping  module  in 
SWAN.  A  slurry  of  China  clay  was  prepared  in  a  tank  by  mixing  clay  powder  with  fresh  tap 
water  in  which  0.5%  sodium  chloride  was  dissolved,  yielding  a  sediment  concentration  of 
about  275  kg/m’.  After  a  consolidation  time  of  six  days  the  concentration  increased  to 
about  500  kg/m3,  i.e.  to  a  density  of  about  1300  kg/m3.  The  viscosity  of  this  suspension 
measured  about  2.7- 1 0'3  m2/s,  i.e.  about  300  times  larger  than  that  of  water.  Though,  the 
strength  of  the  mud  has  not  been  measured,  it  is  estimated  at  about  the  liquid  limit. 

Prior  to  test  3,  two  other  tests  were  carried  out.  During  test  1 ,  only  waves  were  generated 
with  a  height  of  22  mm,  whereas  during  test  2  also  a  net  flow  occurred  (velocities  of  5,  10 
and  1 5  cm/s)  at  a  wave  height  of  38  mm.  Then  the  bed  was  left  at  rest  for  another  day  upon 
which  test  3  was  carried  out  with  waves  of  45  mm,  but  no  current. 

Fig.  5.1  shows  the  measured  variation  of  wave  height  at  six  locations  along  the  flume  after 
equilibrium  was  attained.  It  is  shown  that  the  wave  height  at  the  end  of  the  test  section  has 
decreased  to  about  35  mm,  i.e.  a  reduction  by  about  20%. 

Next,  the  new  SWAN  mud  wave  damping  model  was  used  to  simulate  these  measurements. 
The  model  parameters  were  obtained  from  the  measurements  by  De  Wit  (1995);  De  Wit 
and  Kranenburg  (1997)  and  are  summarized  in  Table  5.1. 
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Fig.  5.1.  Comparison  of  measured  and  computed  wave  height  in  laboratory  flume. 


The  results  of  the  SWAN-simulations  are  also  presented  in  Fig.  5.1,  showing  a  favorable 
agreement  between  the  experimental  data  and  the  numerical  simulations. 


Table  5.1 


A  [kg/m3! 

Pm  fkg/nr  1 

Sm.o  fm] 

h0  [m] 

H[m] 

vm  [nr/sl  a 

1000 

1300 

0.115 

0.325 

0.045 

2-6. 1 0*3  1 

5.2  Application  to  Guyana  coast 

5.2. 1  Observations 

Fig.  5.2  shows  an  aerial  photograph  of  the  breaker  zone  at  Demerara  coast  of  the  Guyana 
coastal  system.  Towards  the  west  (right-hand  part  of  the  photo)  white  bands  of  foam, 
characteristic  for  wave  breaking,  are  clearly  visible.  Further  to  the  east,  no  foam  is  visible; 
in  fact  hardly  any  waves  can  be  observed. 

However,  no  detailed  wave  data  in  the  Guyana  coastal  system  are  available  at  present  from 
which  the  wave  attenuation  over  the  muddy  seabed  can  be  assessed  quantitatively.  Only 
some  general  data  have  been  reported  by  NEDECO  (1972)  on  decreases  in  wave  heights 
under  average  conditions  from  1 .3  to  0.3  m  and  under  extreme  conditions  from  4  to  1 .2  m. 


Fig.  5.2:  Aerial  photograph  of  wave  braking  and  ware  damping  at  Guyana,  Demerara  coast. 


A  good  qualitative  picture  can  be  obtained  from  data  a  similar  site  though.  Fig.  5.3  depicts 
the  measurements  by  Wells  and  Coleman  (1986)  of  wave  spectra  at  three  locations  in  the 
Surinam  coastal  waters.  The  spectra  consist  of  two  parts:  a  strong  peak  at  Tw  =  9  -  10  s, 
where  T„  is  the  wave  period,  reflecting  energy  in  the  swell  from  the  Atlantic  Ocean  and  a 
wider  band  between  Tw  =  ~2  and  -5  s,  reflecting  the  energy  in  locally  generated  waves.  The 
total  wave  energy  has  dissipated  by  about  88%  at  the  middle  measuring  station  (i.e.  over  a 
distance  of  about  9  km)  and  for  96%  at  the  most  near-shore  station.  Note  the  difference  in 
vertical  scale  in  the  graphs.  The  conditions  at  Surinam  coast  are  so  similar  to  those  at  the 
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Guyana  coastal,  that  it  is  believed  that  the  observations  of  Fig.  5.3  are  characteristic  for  the 
Guyana  coastal  system  as  well. 

As  explained  before,  the  strong  wave  attenuation  shown  in  Fig.  5.2  and  5.3  is  expected  to 
be  caused  by  viscous  dissipation  in  the  liquefied  seabed.  Unfortunately,  up  till  now  we 
have  no  direct  measurements  of  such  liquefaction,  and  we  have  to  rely  on  physical 
interpretation  of  other  data.  Fig.  5.6  shows  the  filtered  results  of  recent  dual  frequency  echo 
soundings  along  a  cross  section  perpendicular  to  Guyana  coast.  The  upper  signal  (210  kHz) 
reflects  the  top  of  the  soft  mud  bed,  whereas  the  lower  signal  (33  kHz)  reflects  some  of  the 
structure  within  the  bed.  This  graph  reveals  a  small  pool  of  fluid  mud  only,  depicted  by  the 
gap  between  the  two  signals  at  around  x  =  8000  m.  The  irregularity  of  the  lower  frequency 
echo  signal  is  typical  for  layered  mud  deposits.  This  is  further  depicted  in  Fig.  5.5,  showing 
reflections  from  a  20  kHz  sub-bottom  profiler  applied  by  Delft  Hydraulics  Laboratory 
(1962)  in  the  vicinity  of  Demerara  river  mouth.  The  multiple  reflections  are  typical  for 
stratified  mud  deposits;  the  thickness  of  the  layers  is  estimated  at  about  1  -  2  dm  (the 
vertical  resolution  is  not  known  unfortunately,  but  can  be  estimated  from  the  amplitude  of 
the  wavy  acoustic  signal3').  This  layered  structure  cannot  be  caused  by  sedimentation 
events,  as  the  sedimentation  rate  in  the  Guyana  coastal  system  is  not  driven  by  episodic 
events.  Hence,  this  stratification  can  only  be  caused  by  frequent,  incomplete  liquefaction  of 
the  seabed. 


31  These  small  acoustic  waves  are  marked  by  an  ellipse  in  Fig.  5.5.  There  amplitude  a  follows  from  a 
=  c/f  where  c  is  the  speed  of  sound  in  the  seabed  (c  «  1500  m/s)  and/is  the  frequency  of  the 
acoustic  wave  (f=  20  kHz).  Hence,  a  ~  7.5  cm. 
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From  this  analysis  we  deduce  from  Fig.  5.4  and  5.5  that  the  thickness  of  the  liquefiable  bed 
amounts  to  about  2  m,  and  we  infer  that  we  can  assess  the  thickness  of  the  fluid  mud  layer 
and  its  extension  from  spatial  data  on  the  echo-soundings. 

The  in-situ  density  has  been  measured  by  NEDECO  (1972).  They  found  values  for  pm 
ranging  from  1300  to  1500  kg/m’.  As  variations  in  density  cannot  be  very  large,  and  as  the 
dissipation  function  is  not  very  sensitive  to  small  variations  in  density,  we  take  an  average 
value  of  pm  =  1400  kg/m3. 

Some  viscosity  measurements  have  been  carried  out  by  NEDECO  (1972),  as  well,  with 
ranges  between  rj  =  0.1  -  10  Pa  s.  As  these  data  are  not  very  accurate,  the  viscosity  is  our 
main  calibration  parameter. 
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Fig.  5.4:  Filtered  acoustic  reflections  from  dual  frequency  echo  sounding  (f  =  210  &  33  kHz)  along  a  cross 
_ section  in  Guyana  coastal  system  perpendiculars  to  the  coast. _ 


X:  337786 
Y:  806774 
0 

20 
10 
18 
14 
12 
10 
8 
6 
4 
2 
0 

Base  Scrvey 
Line  2033  . 


Line  2033  .. 

2000  4000 


Length  14866  Azimuth:  90 

6000  8000  10000 


X-  352652 
V  806790 

12000  14000 


Fig.  5.5:  Acoustic  reflections  from  sub-bottom  profiling  (f  ~  20  kHz)  in  Guyana  coastal  system,  Demerara 
mouth  ( Delft  Hydraulics  Laboratory,  1962)  -  the  vertical  scale  is  not  know  n,  but  estimated  at  about  0.5  m. 


5.2.2  Simulations 

The  modified  SWAN-model  (as  described  in  Chapter  3)  is  applied  to  the  wave  model 
covering  a  large  part  of  the  Guyana  coastal  system,  including  the  Essequibo  estuary  (see 
Fig.  5.6).  The  model  covers  a  250  km  wide  coastal  area  and  extends  about  85  km  offshore. 
The  depth  at  the  offshore  boundary  of  the  model  varies  between  135  to  40  m. 
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Fig.  5.6:  Depth  scheniatisation  of  Guyana  wave  model. 


The  wave  model  employs  a  curvilinear  grid  that  comprises  263  by  383  points,  with 
approximately  42,000  active  points.  This  computational  grid  is  presented  in  Fig.  5.7.  The 
dimensions  of  the  grid  cells  vary  from  2.5  km  near  the  model  boundaries  to  approximately 
200  m  in  the  Essequibo  estuary,  near  Georgetown  and  in  the  nearshore  zone. 
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Fig.  5.7.  Computational  curvilinear  grid  of  Guyana  Wave  model. 
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Considering  the  presence  of  relatively  long  waves,  the  width  of  the  energy  distribution  was 
taken  at  ms  =  10  (instead  of  ms  =  4  in  case  of  wind-generated  waves)  in  the  wave  model. 
For  the  wave  spectrum  a  directional  resolution  of  10°  was  adopted  over  the  full  circle.  The 
frequency  spreading  was  schematised  using  3 1  discrete  frequencies  between  0.05  Hz  and  1 
Hz  using  a  logarithmic  scale.  For  the  remaining  parameters  the  default  settings  of  the  SWAN 
model  were  used. 

The  swan  mud  wave  damping  model  requires  as  input  the  thickness  of  the  (fluid)  mud 
layer  in  the  model  area.  This  was  defined  as  the  difference  between  the  high  frequency 
echo  reflections  and  the  lowest  10%  of  the  low  frequency  echo  reflections  (see  Fig.  5.4). 
The  resulting  spatial  varying  mud  layer  thickness  varying  between  1  m  and  5.5  m,  as 
indicated  in  Fig.  5.8. 
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Fig.  5.8.  Mud  lay  er  thickness  as  applied  in  SWAN  mud  wave  damping  model. 
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The  model  is  calibrated  on  the  basis  of  a  comparison  of  average  wave  heights  observed  at 
three  nearshore  locations,  as  presented  in  NEDECO  (1972).  The  observation  locations  are 
indicated  in  Figure  5.6.  The  objective  of  the  calibration  is  to  optimise  the  settings  of  the 
input  variables  of  the  SWAN  mud  wave  damping  model.  A  further  calibration  of  the  model 
on  the  basis  of  simultaneous  offshore  and  nearshore  measurements  is  not  possible  as  no 
simultaneous  nearshore  wave  measurement  data  are  available  at  present. 
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In  NEDECO  (1972)  average  and  maximum  wave  heights  at  three  nearshore  and  one  offshore 
locations  are  presented.  Considering  the  dynamic  behaviour  of  the  (mud)  bed  at  the  location 
of  the  nearshore  observations,  it  is  very  likely  that  the  current  seabed  differs  from  the  seabed 
at  the  time  the  observations  were  made.  Consequently,  it  is  not  possible  to  make  a  direct 
comparison  between  the  nearshore  observations  and  model  results  (based  on  the  current 
bathymetry).  It  is  therefore  not  possible  to  accurately  represent  these  observations. 
However,  we  require  them  to  be  of  the  same  order  of  magnitude  as  an  indication  to  whether 
the  damping  of  the  wave  energy  is  properly  represented  by  the  SWAN  mud  wave  damping 
model. 

As  no  detailed  information  was  available  on  the  characteristics  of  the  fluid  mud  itself,  the 
mud  parameters  had  to  be  assessed  from  literature  or  sensitivity  analyses.  We  kept  the 
density  of  the  fluid  mud  constant,  as  this  is  expected  not  to  vary  too  much.  The  mud’s 
viscosity  is  determined  by  trial-and-error  and  is  thus  basically  a  calibration  factor.  The 
selected  values  are  given  in  Table  5.2.  Computations  were  carried  out  with  vm  =  0.01,  0.001 
and  0.0001  m2/s. 

Table  5.2 


Input  parameters  and  settings  for  the  SWAN  mud  wave  damping  model. 


Parameter 

settings 

thickness  of  the  fluid  mud  layer  Sw 

spatially  varying  on  computational  grid  [m] 

density  of  fluid  mud,  pft1 

1400  kg/m3 

density  of  water,  pw 

1010  kg/m'  (at  T  =  25°C  and  5=  20  ppt) 

viscosity  of  fluid  mud,  Vm 

calibration  parameter  [nf/s],  ranging  between  0.01  and  0.0001  nf/s 

calibration  parameter,  a 

1.0  [-] 

The  effect  of  the  varying  v,„  during  an  extreme  wave  condition  (H„,o  =  4  m,  Tp  =  10  s 
accompanied  by  a  NE  wind  of  12.5  m/s)  on  the  significant  wave  height  is  presented  in  Fig. 
5.9.  This  figure  shows  cross-shore  profiles  (approximately  located  at  Hope  Beach,  see  Fig. 
5.6)  of  the  significant  wave  height  and  depth  profile  as  well  as  the  thickness  of  the  fluid  mud 
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layer  as  schematised  in  the  SWAN  model.  It  shows  that  the  mud’s  viscosity  has  an  important 
effect  on  the  wave  attenuation. 


In  our  further  analyses  we  apply  vm  =  10'1  nr/s  (i.e.  rj„,  =  1.3  Pa  s),  as  this  value  gave  the 
best  results;  moreover,  this  is  a  typical  value  for  muddy  beds  at  the  observed  densities.  The 
effect  of  fluid  mud-induced  wave  damping  in  the  wave  computations  is  shown  in  Figure 
5.10  for  an  extreme  wave  conditions  from  the  ENE  direction.  This  figure  presents  cross¬ 
shore  profiles  of  the  significant  wave  height  peak  and  mean  wave  period  (Tp  and 
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T„,oi );  wave  dissipation  m  )  ;  and  depth  profile  and  thickness  of  the  fluid  mud  layer  (ho 
and  Sm  o)  as  schematised  in  the  model.  Figure  5.9  shows  that: 

The  mean  wave  period  decreases  strongly  over  the  fluid  mud  layers,  with  a  stronger 
dissipation  of  the  low-frequency  waves; 

The  wave  heights  start  to  decrease  a  about  25  m  water  depth  with  a  relative  decrease  up  to 
50%  at  about  5  m  water  depth  when  taking  into  account  wave  damping  by  fluid  mud; 

The  dissipation  of  wave  energy  starts  at  the  offshore  end  of  the  fluid  mud  layer  and  is  fairly 
constant  towards  the  coast.  In  case  wave  damping  by  fluid  mud  is  not  included  in  the  wave 
model,  a  peak  in  wave  dissipation  is  found  at  the  coastline,  where  the  waves  break  due  to 
depth  limitation.  The  latter  not  being  observed  at  the  area  of  interest  of  Guyana  coastal 
system. 

The  effect  of  wind-generated  waves  on  the  wave  climate  is  illustrated  in  Fig.  5.1 1  showing 
one-dimensional  wave  spectra  at  three  locations  for  an  extreme  wave  condition  (Hs  =  4  m) 
from  the  NE  direction.  The  simulation  includes  wave  damping  by  fluid  mud.  Note  the 
growth  of  wind-generated  waves  (at  a  frequency  of  about  0.2  Hz)  towards  the  coast  and  the 
dissipation  of  the  swell  (at  a  frequency  of  about  0. 1  Hz).  In  this  case  the  peak  period  shifts 
from  10  s  offshore  to  5  s  at  a  nearshore  location.  The  energy  of  the  >  0.2  Hz  waves  reduces 
by  a  factor  3  from  A  to  B,  and  at  their  peak  by  a  factor  5.10.  This  uneven  behaviour  in 
dissipation  is  owing  to  the  effects  of  bed  friction  and/or  viscous  dissipation  in  the  mud  layer. 
The  latter  effect,  however,  is  by  far  the  larger  of  the  two  (e.g.  Fig.  5.10).  This  process  of 
wave  attenuation  towards  the  coast  agrees  qualitatively  with  the  measurements  by  Wells  and 
Coleman  (1986),  see  Fig.  5.3. 


Fig.  5.1 1:  Computed  wave  spectrum  at  three  locations  in  Demarara  coastal  system  with  ( black  line )  and 
_ without  (red  line)  locally  generated  waves  (note  different  scales). _ 
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Table  5.3  presents  the  comparison  between  historical  nearshore  observations  described  in 
NEDECO  (1972)  and  the  model  results  with  and  without  the  fluid  mud-induced  wave 
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damping.  Here  it  is  concluded  that  when  including  the  swan  mud  wave  damping  model 
(with  the  above  variable  settings),  the  damping  effect  of  the  wave  energy  is  in  reasonable 
agreement  with  the  observations.  Wave  heights  are  significantly  being  overestimated  in  case 
fluid  mud-induced  wave  damping  is  not  included  in  the  SWAN  simulations. 

Table  5.3 

Comparison  between  historical  observations  (NEDECO,  1972)  and  model  results  with  (Dm 
=  0.001  m2/s. 


Location 

observed  Hs  [m] 

Computed  Hs  [m] 

Computed  Hs  [m] 

with  fluid  mud 

without  fluid  mud 

Average  wave  conditions 

Offshore 

1.3 

1.3 

1.3 

Buxton 

0.6 

0.6 

1.0 

Demerara 

0.5 

0.6 

1.0 

Beacon 

Kitty  Groyne 

0.3 

0.1 

0.5 

Maximum  wave  conditions 

Offshore 

4.0 

4.0 

4.0 

Buxton 

2.0 

1.6 

2.9 

Demerara 

1.3 

1.5 

2.8 

Beacon 

Kitty  Groyne 

1.2 

0.4 

0.9 
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6  Application  of  wave-induced  mud  transport 


Sakakiyama  and  Bijker  (1987)  describe  experiments  in  a  horizontal  wave  flume,  measuring 
both  mud-induced  wave  damping  and  wave-induced  mud  transport.  A  bed  of  kaolinite  and 
fresh  water  was  placed  in  a  test  section  in  between  two  wooden  false  floors.  The  test 
section  was  12  m  long  and  the  mud  thickness  amounted  0.095  m,  whereas  the  flume  was 
0.5  m  wide.  Mud  beds  with  a  density  varying  between  1 140  and  1380  kg/m’  were  prepared. 
Wave  heights  were  varied  between  1  and  4.4  cm,  at  wave  periods  varying  between  0.6  and 
2  s.  The  water  depth  was  kept  constant  at  0.3  m. 

Sakakiyama  and  Bijker  (1987)  measured  surface  wave  height  and  attenuation,  the  height  of 
the  water-mud  interface,  and  the  velocity  within  the  mud  layer.  Also,  viscosity  and  strength 
of  the  mud  layer  were  established  independently. 

The  series  at  the  highest  mud  density  are  most  interesting  as  they  show  considerable  wave 
damping  at  higher  wave  heights  and  ditto  mass  transport  in  the  first  part  of  the  mud  bed, 
whereas  in  the  second  part  of  the  mud  bed,  waves  are  damped  so  much  that  liquefaction 
nor  mass  transport  was  measured.  This  is  depicted  in  Fig.  6.1,  showing  for  the  higher  initial 
wave  height  damping  over  the  first  7  m  of  the  test  section,  where  after  the  wave  height 
remains  constant  at  about  1  cm.  When  the  initial  wave  height  amounted  to  1  cm,  no 
damping  occurred;  apparently,  the  mud  bed  could  not  be  liquefied. 


Fig.  6.  / :  Wave  attenuation  measured^  by  Sakakiyama  and_  Bijker  II 9891 


When  mobile,  the  velocity  within  the  mud  bed  amounted  to  about  0.1-1  mm/s,  which  was 
measured  at  about  x  =  3.5  m.  The  vertical  velocity  profiles  were  very  homogeneous  with 
little  vertical  gradients. 

Fig.  6.2  shows  very  little  variation  of  mud  velocity  with  wave  height.  This  implies  a 
balance  between  the  wave-induced  traction  and  friction  (yield  stress),  where  the  traction  is 
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a  linear  function  of  the  dissipation  rate  and  gradient  in  wave  energy.  The  dissipation  rate 
appears  to  be  more  or  less  constant  along  the  flume  over  the  section  where  waves  are 
damped,  e.g.  Fig.  6.1. 


wave  height  H  [cm] 


Fig.  6. 2:  Fluid  mud,  velocities  measure d  by  Sakakivama  amj_  Bicker  [ /  989). 


The  parameters  describing  Sakakiyama  and  Bijker’s  series  A  experiment  are  given  in  Table 

6.1. 


Table  6.1 


A  [kg/nr] 

Pm  fkg/nr  1  Sm.„  [ml 

ho  [m] 

H[  m] 

T[s] 

V„  [nr/sl  a 

1000 

1370  0.09 

0.3 

0.035 

1.0 

(1  -  2)- 1  O'2  1 

Dus  test  sow 

bestaat  uit  eerste  5  m 

van  goot ,  conform 

Fig.  6.1. 

Narekenen  van  het  hele 

experiment  is  zinnig  als  we  ook  rheologisch  model  hebben  ontwikkeld. 


Eerste  serie  numerieke  experimenten: 

The  first  series  of  numerical  experiments  to  test  the  combined  effects  of  fluid  mud  induced 
wave  damping  and  wave-induced  fluid  mud  transport  is  carried  out  on  a  horizontal  or 
sloping  bed  with  waves  at  perpendicular  direction  -  this  simulates  a  coast  with  cross-shore 
waves.  This  lay-out  is  sketched  in  Fig.  6.3. 

We  carry  out  two  series  of  experiments.  Series  1:  now  waves,  series  2:  with  waves.  In  both 
series,  the  slope  of  the  bed  i  is  increased  in  small  steps  from  0  to  1  O  '  rad.,  and  in  both 
series  there  is  no  water  flow  in  the  upper  layer.  Moreover,  erosion,  entrainment, 
sedimentation  and  consolidation  are  not  included.  The  entire  length  of  the  model  amounts 
to  2.5  km,  with  20  m  horizontal  resolution.  The  width  of  the  model  is  set  to  2  km,  and  the 
width  of  the  fluid  mud  layer  at  1  km.  Table  6.2  gives  an  overview  of  the  various 
parameters: _ 
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Table  6.2 


mud  parameters  in  2L  FLUID  MUD  model  for  Series  1  &  2  experiments. 


parameter 

value 

comments 

ho 

5  m 

largest  depth  in  model  -  shallower  “towards  beach" 

L 

2500  m 

B 

15,000  m 

very  wide  because  of  wave  spreading 

5m 

0.2  m 

Lm 

60  (60)  600  m 

to  be  varied 

Ls 

>  150  m 

Bm 

600  m 

symmetrical  within  overall  domain 

Ax  =  Ay 

20  m 

At 

2  min 

based  on  internal  celerity  and  Cr  =  5; 
to  be  varied  (A/  =  0.5;  1 ;  5;  10  min) 

i 

0  (0.2- 10'3)  1 0'3  rad 

to  be  varied 

Hs 

2.0  m 

T 

8  s 

A 

1000  kg/m3 

Pm 

1370  kg/m3 

Vm 

1  - 1 0 2  m2/s 

tb 

1  Pa 

A 

0.02 

corresponds  to  C~  60  m'Vs 
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The  time  step  follows  form  the  internal  wave  speed  c,  =  yj(pm~  PK ) £g8m  / pw  and  a 

Courant  number  not  larger  than  about  5.  Hence,  At  <  5 Ax/Cj,  i.e.  1 17  s  for  the  values  of 
Table  6.2.  Thus  At  is  set  to  2  minutes  as  a  reference. 


The  configuration  of  Fig.  6.3  needs  some  refinement,  as  the  density  gradient  at  the  edges  of 
the  fluid  mud  is  too  large  to  remain  stable  conditions  in  case  of  even  a  horizontal  bed. 
Hence,  edges  as  sketched  in  Fig.  6.3  should  be  added.  Their  length  can  be  deduced  from 
the  balance  between  density  gradient  and  bed  shear  stress  (equ.  (3.5)): 


gAp  8Zm  TZmJ  =Q 

Pm  dx,  PJ 


(6.1) 


This  leads  to  Ls  >  gd'Ap/ tB  ,  amounting  to  147  m  for  the  values  of  Table  6.2. 


The  wave  length  for  these  conditions  amounts  to  about  32  m,  hence  the  real  wave  number 
k,  =  0.2  m'1  and  kho  =  0.98.  We  have  chosen  the  conditions  such  that  the  imaginary  wave 
number  is  near  its  maximum,  estimated  at  k ,  ~  0.001  m'1.  Hence,  within  150  m,  95%  of  the 
wave  energy  is  expected  to  be  dissipated  within  the  fluid  mud  layer. 

The  Bingham  strength  Tb  is  chosen  such  that  the  fluid  mud  should  start  to  slide  downhill  at 
a  bed  slope  i  >  0.5  10‘?  rad. 

What  effect  of  waves  can  be  expected?  For  waves  perpendicular  to  the  coast  (in  x- 
direction),  6  -  0,  and  the  radiation  stress  in  this  direction  for  monochromatic  waves 
amounts  to 


T,=-^(pg{2,,-\l2)E)P-^(pg{n-V2)H;)*f^  {6, 

As  2/7  -  0.5  =  0{  1}.  We  have  selected  the  wave-mud  conditions  such  that  S,„{ d2  v)05  = 
0.2(71/8/0.01)" 5  =  1.25,  which  yields  about  maximum  damping.  As 

Hs  =  Hs  0  expj-^jr} dH2 /dx  =  -2ktH2H  Substitution  in  (6.2)  and  use  of  (3.42)  yields: 


P  8 

r  m  m 


*  7 


(6.3) 


From  our  SWAN  simulations  we  found  that  kj  ~  0.5- 10°  m  '.  Further,  r)2 /t)2  ~  0.1,  and  we 
Tr,m,x / Pm8m  *  0-6  •  1 0"  H2  (=  0.0024  for  Hs  =  2  m).  This  acceleration  is  small 

compared  to  gravitational  effects  of  the  mud  on  a  slope.  For  instance  a  =  0.001  one  finds 
gApa/ p  =  0.0036 .  Only  for  waves  with  Hs  >  2.5  m,  the  radiation  stresses  will  exceed 
the  gravitational  effects. 

However,  in  that  case,  the  Bingham  strength  is  still  working,  and  the  difference  between 
the  radiation  stresses  and  gravitational  effects  for  the  conditions  under  consideration  will 
never  be  large  enough  to  overcome  the  Bingham  strength  and  move  the  fluid  mud  uphill. 
This  implies  that  this  is  only  possible  when  the  mud  is  truly  liquefied  and  has  lost  all  its 
strength. 
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7  Application  of  wave-induced  liquefaction 


To  be  completed 


7-  I 


WL  |  Delft  Hydraulics 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud-induced  wave  transport  &  wave-induced  mud 
transport 


Implementation  in  DELFT3D 


8  References 


Alba,  P.D.,  Seed,  H.B.,  Chan,  C.K.,  1976,  Sand  liquefaction  in  large-scale  simple  shear  tests,  ASCE, 
Journal  of  Geotechnical  Engineering  Division,  Vol  102  (GT9)  909-927. 

Allison,  M.A.  and  M.T.  Lee,  2004,  Sediment  exchange  between  the  Amazon  mud  banks  and  shore- 
fringing  mangroves  in  French  Guiana,  Marine  Geology,  208,  193-203. 

Augustinus,  P.G.E.F.,  204,  The  influence  of  the  Trade  Winds  on  the  coastal  development  of  the 
Guianas  a  various  scale  levels:  a  synthesis.  Marine  Geology,  208,  145-153. 

Booij,  N.,  R.C.  Ris  and  L.FI.  Holthuijsen,  1999,  A  third-generation  wave  model  for  coastal  regions. 
Part  I:  Model  description  and  validation.  Journal  of  Geophysical  Research,  104  (C4)  7649-7666. 

Chou,  H.-T.,  1989,  Rheological  response  of  cohesive  sediments  to  water  waves,  PhD-dissertation, 
University  of  California,  Berkeley,  USA. 

Chou,  H.-T.,  M.A.  Foda  and  J.R.  Hunt,  1993,  Rheological  response  of  cohesive  sediments  to 
oscillatory  forcing,  in:  Coastal  and  Estuarine  Studies,  Vol  42,  Nearshore  and  Estuarine 
Cohesive  Sediment  Transport,  ed.  A.J.  Mehta,  126-147. 

Coussot,  P.,  1997,  Mudflow  rheology  and  dynamics,  IAHR  monograph  series,  Balkema,  Rotterdam, 
pp  255. 

Dalrymple,  R.A.  and  P.L.  Liu,  1978,  Waves  over  soft  mud  beds:  a  two-layer  fluid  mud  model. 
Journal  of  Physical  Oceanography,  8,  1121-1131. 

Delft  Hydraulics  Laboratory,  1962,  Demerara  coastal  investigation;  Report  on  siltation  of  Demerara 
Bar  Channel  and  coastal  erosion  in  British  Guyana,  Vol  2,  Field  survey  and  basic  studies.  The 
Netherlands. 

Delft  Hydraulics  Laboratory,  1974,  Momentum  and  mass  transfer  in  stratified  flows.  Report  R880. 

De  Wit,  P.J.,  1995,  Liquefaction  of  cohesive  sediment  by  waves,  PhD-dissertation,  Delft  University 
of  Technology,  The  Netherlands. 

De  Wit,  P.J.  and  C.  Kranenburg,  1997,  On  the  liquefaction  and  erosion  of  mud  due  to  waves  and 
current,  in:  Cohesive  Sediments,  ed.  N.  Burt,  R.  Parker  and  J.  Watts,  John  Wiley  &  Sons,  33 1  - 
340. 

Dingemans,  M.W.,  1997,  Water  wave  propagation  over  uneven  bottoms;  Part  1  &  II,  World 
Scientific,  Singapore. 

Foda,  M.A.  and  S.-Y.  Tzang,  1994,  Resonant  waves  of  silty  soil  by  water  waves.  Journal  of 
Geophysical  Research,  99  (CIO)  20,463-20,475. 

Foda,  M.A.,  J.R.  Hunt  and  H.-T.  Chou,  1993,  A  nonlinear  model  for  the  fluidization  of  marine  mud 
by  waves,  Journal  of  Geophysical  Research,  16  (2)  7039-7047. 

Gade,  H.G.,  1958,  Effects  of  a  non-rigid,  impermeable  bottom  on  plane  surface  waves  in  shallow 
water.  Journal  of  Marine  Research,  16  (2)  61-82. 

Golden,  S.P.,  J.W.  Godwin  and  A.D.  Olal,  1982,  The  dependence  of  the  elastic  properties  of  clay 
dispersions  on  the  mode  of  interaction  between  the  particles.  Transactions  and  Journal  of  the 
British  Ceramic  Society,  81  (3)  84-87. 

Groeneweg,  J.,  1999,  Wave-current  interactions  in  a  Lagrangean  mean  formulation,  PhD-thesis, 
Delft  University  of  Technology,  Department  of  Civil  Engineering  and  Geosciences. 

Healy,  T.,  Y.  Wang  and  J.-A.  Healy,  2002,  Muddy  coasts  of  the  world:  Processes,  Deposits  and 
Function,  Elsevier,  Proceedings  in  Marine  Science  No  4,  542  pp. 

Jiang,  L.  and  Z.  Zhao,  1989,  Viscous  damping  of  solitary  waves  over  fluid-mud  beds,  ASCE, 
Journal  of  Waterways,  Port,  Coastal  and  Ocean  Engineering,  1 15  (3)  345-362. 

Jiang,  F.,  1993,  Bottom  mud  mass  transport  due  to  water  waves,  PhD-thesis,  University  of  Florida, 
Gainesville,  Florida,  USA. 

Kump,  1972,  Littoral  currents  in  relation  to  the  mud  bank  formation  along  the  coast  of  Kerala, 
Mahasagar,  Bulletin  of  the  National  Institute  of  Oceanography,  5  (3). 


8-  I 


WL  |  Delft  Hydraulics 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud-induced  wave  transport  &  wave-induced  mud 
transport 


Implementation  in  DELFT3D 


Lambe,  T.W.  and  R.V.  Whitman,  1979,  Soil  Mechanics,  SI  Version,  John  Wiley  and  Sons,  New 
York,  pp  553. 

Lindenberg,  J.,  L.C.  Van  Rijn  and  J.C.  Winterwerp,  Some  experiments  on  wave-induced 
liquefaction  of  soft  cohesive  soils,  Journal  of  Coastal  Research,  5,  127-138. 

Liu,  P.L.-F.,  1973,  Damping  of  water  waves  over  porous  bed,  ASCE,  Journal  of  the  Hydraulic 
Division,  92  (12)  2263-2271. 

Liu,  K.  and  C.C.  Mei,  1989,  Effects  of  wave-induced  friction  on  a  muddy  seabed  modelled  as  a 
Bingham-plastic  fluid.  Journal  of  Coastal  Research,  5  (4)  777-789. 

Longuet-Higgins,  M.S.,  1953,  Mass  transport  in  water  waves.  Philosophical  Transactions  of  the 
Royal  Society  of  London,  A  (903)  535-581. 

Longuet-Higgins,  M.S.,  1958,  The  mechanics  of  the  boundary-layer  near  the  bottom  in  a  progressive 
wave,  in:  Proceedings  of  the  6th  Conference  on  Coastal  Engineering,  Gainesville,  Florida,  ed. 
J.W.  Johnson,  Council  on  Wave  Research;  The  Engineering  Foundation,  184-193. 

Maa,  P.-Y.,  1986,  Erosion  of  soft  mud  beds  by  waves,  PhD-dissertation,  University  of  Florida, 
Coastal  and  Oceanographic  Engineering  Department,  Gainesville,  Florida,  USA,  rep. 
UFL/COEL-TR-059. 

Maa,  P.-Y.  and  A.J.  Mehta,  1990,  Soft  mud  response  to  water  waves,  ASCE,  Journal  of  Waterways, 
Port,  Coastal  and  Ocean  Engineering,  1 16  (5)  634-650. 

MacPherson,  H.,  1980,  The  attenuation  of  water  wave  over  a  non-rigid  bed.  Journal  of  Fluid 
Mechanics,  97  (4)  721-742. 

Mallard,  W.W.,  Dalrymple,  R.A.,  1977,  Water  waves  propagating  over  a  deformable  bottom. 
Proceedings  of  the  9th  Offshore  Technology  Conference,  Houston,  TX,  Vol  3,  141-146. 

Malvern,  L.E.,  1969,  Introduction  to  the  Mechanics  of  a  Continuous  Medium,  Prentice-Hall,  Series 
in  Engineering  of  the  Physical  Sciences,  New  Jersey. 

McDougal,  W.G.,  Y.T.  Tsai,  P.-F.  Liu,  E.C.  Clukey,  1989,  Wave-induced  pore  water  pressure 
accumulation  in  marine  soils.  Journal  Offshore  Mechanics  and  Artie  Engineering,  ASME,  111, 
1-11. 

Mei,  C.C.  and  K.F.  Liu,  1987,  A  Bingham  plastic  model  for  a  muddy  seabed  under  long  waves. 
Journal  of  Geophysical  Research,  94  (Cl  3)  14,581-14,594. 

Mei,  C.C.,  1989,  Applied  dynamics  of  ocean  surface  waves.  World  Scientific. 

Mehta,  A.J.  1996,  Interaction  between  fluid  mud  and  water  waves,  in:  Environmental  Hydraulics, 
ed.  V.P.  Singh  and  W.H.  Hager,  Kluwer  Academic  Publishers,  The  Netherlands,  153-187. 

Moore,  F.,  1959,  The  rheology  of  ceramic  slips  and  bodies.  Transactions  of  the  British  Ceramic 
Society,  58,  470-494. 

NEDECO,  1972,  Report  on  sea  defense  studies  of  British  Guyana,  NEDECO,  The  Hague,  The 
Netherlands. 

Phillips,  O.M.,  1977,  The  dynamics  of  the  upper  ocean,  Cambridge  University  Press. 

Rodriguez,  H.N.,  2000,  Mud  bottom  evolution  at  open  coasts,  PhD-dissertation,  University  of 
Florida,  Coastal  and  Oceanographic  Engineering  Department,  Florida,  USA. 

Sakakiyama,  T.,  E.W.  Bijker,  1989,  Mass  transport  velocity  in  mud  layer  due  to  progressive  water 
waves,  ASCE,  Journal  of  the  Waterway,  Port,  Coastal,  and  Ocean  Engineering,  115  (5)  614- 
633. 

Sheremet,  A.,  G.W.  Stone,  2003,  Observations  of  nearshore  wave  dissipation  over  muddy  sea  beds. 
Journal  of  Geophysical  Research,  1 08  (C 1 1 ),  2 1  - 1  -2 1  - 1 1 . 

Soulsby,  R.L.,  S.  Clarke,  2005,  Bed  shear-stresses  under  combined  waves  and  currents  on  smooth 
and  rough  beds,  HR  Wallingford,  Report  TR  137. 

Spierenburg,  S.E.J.,  1987,  Seabed  response  to  water  waves,  PhD-dissertation,  Delft  University  of 
Technology,  The  Netherlands. 

Stelling  G.S.,  Kernkamp,  H.W.J.  and  Laguzzi,  M.M.,  1998,  ‘‘Delft  Flooding  System:  a  powerful  tool 
for  inundation  assessment  based  upon  a  positive  flow  simulation",  in  Hydro-informatics'98,  ed. 
by  Babovic  &  Larsen,  Balkema,  Rotterdam,  pp  449-460. 


WL  |  Delft  Hydraulics 


8-2 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud-induced  wave  transport  &  wave-induced  mud 
transport 


Implementation  in  DELFT3D 


Sumer,  B.M.,  J.  Fredsoe,  2002,  The  mechanics  of  scour  in  the  marine  environment.  World 
Scientific,  Advanced  Series  on  Ocean  Engineering,  Vol.  17. 

Toorman,  E.A.,  1997,  Modelling  the  thixotropic  behaviour  of  dense  cohesive  sediment  suspensions, 
Rheologica  Acta,  36,  56-65. 

Tsuruya,  H.,  S.  Nakano,  J.  Takahama,  1987,  Interaction  between  surface  waves  and  a  multi-layered 
mud  bed.  Report  of  the  Port  and  Harbour  Research  Institute,  Japan,  26  (5)  137-173. 

Tubman,  M.W.,  Suhayda,  J.N.,  1976,  Wave  action  and  bottom  movements  in  fine  sediments,  Proc. 
15th  ASCE  conference,  1168-1183. 

Van  Kesteren,  W.G.M.,  Comelisse,  J.M.,  1995,  Modeling  of  the  behavior  of  natural  mud  beds  under 
cyclical  or  wave  loading,  Rijkswaterstaat  and  Delft  Hydraulics,  Research  on  Cohesive 
Sediment,  Report  50. 

Van  Kesteren,  W.G.M.,  2000,  Slope  stability  analysis  Orange  River  mouth,  WL  |  Delft  Hydraulics, 
Report  Z2739. 

Verbeek,  H.  and  J.M.  Comelisse,  1997,  Erosion  and  liquefaction  of  natural  mud  under  surface 
waves,  in:  Cohesive  Sediments,  ed.  N.  Burt,  R.  Parker  and  J.  Watts,  John  Wiley  &  Sons,  353- 
364. 

Walstra,  D.J.R.,  D.A.  Roelvink  and  J.  Groeneweg,  2000,  Calculation  of  wave-driven  currents  in  3D 
mean  flow  model,  ****** 

Wang,  Z.B.  and  J.C.  Winterwerp,  1992,  Transport  of  Fluid  Mud,  Numerical  modeling  with  a  two- 
layer  system  -  research  documentation.  Report  Z163. 

Wells,  J.T.  and  J.M.  Coleman,  1981,  Physical  processes  and  fine-grained  sediment  dynamics,  coast 
of  Surinam,  South  America,  Journal  of  Sedimentary  Petrology,  1053-1068. 

Wells,  J.T.  and  G.P.  Kemp,  1986,  Interaction  of  surface  waves  and  cohesive  sediments:  field 
observations  and  geologic  significance,  in:  Lecture  Notes  on  Coastal  and  Estuarine  Studies,  Vol 
14,  Estuarine  Cohesive  Sediment  Dynamics,  ed.  A.J.  Mehta,  43-65. 

Winterwerp,  J.C.,  Z.B.  Wang,  J.A.Th.M.  van  Kester,  J.F.  Verweij,  2003,  On  the  far-field  impact  of 
Water  Injection  Dredging,  ICE/IAHR,  Water  and  Maritime  Engineering,  Vol  154,  Issue  4,  pp  285- 
296. 

Winterwerp,  J.C.  and  W.G.M.  Van  Kesteren,  2004,  An  introduction  to  the  physical  processes  of 
cohesive  sediment  in  the  marine  environment,  Elsevier,  Developments  in  Sedimentology,  56. 

Yamamoto,  T.,  H.L.  Koning,  H.  Sellmeijer  and  E.  Van  Hijum,  1978,  On  the  response  of  a  poro- 
elastic  bed  to  water  waves.  Journal  of  Fluid  Mechanics,  87  (1)  193-206. 

Yamamoto,  T.  and  S.  Takahashi,  1985,  Wave  damping  by  soil  motion,  ASCE,  Journal  of  the 
Waterway,  Port,  Coastal  and  Ocean  Engineering,  111  (1)  62-77. 

Zijlema,  M.  and  A.J.  van  der  Westhuijsen,  205,  On  convergence  behaviour  and  numerical  accuracy 
in  stationary  in  stationary  SWAN  simulations  of  nearshore  wind  wave  spectra.  Coastal 
Engineering,  52,  237-265. 


WL  J  Delft  Hydraulics 


8-3 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud-induced  wave  transport  &  wave-induced  mud 
transport 


Implementation  in  DELFT3D 


WL  |  Delft  Hydraulics 


8-  I 


Modeling  of  mud-wave  interaction 


Z3672 


November  2007 


Mud-induced  wave  transport  &  wave-induced  mud 
transport 


Implementation  in  DELFT3D 


A  Appendix:  List  of  symbols 


A  wave  semi-orbital  excursion 

a  aggregation  coefficient  in  rate  of  flocculation  equation 

b  break-up  coefficient  in  rate  of  flocculation  equation 
Cs  coefficient  in  entrainment  model 

CG  coefficient  in  entrainment  model 

c  sediment  concentration  by  mass  in  upper  layer 
c  cohesion 

Co  cohesion  at  fully  flocculated  conditions 
c  wave  celerity 

cg  group  velocity 

cm  sediment  concentration  by  mass  in  fluid  mud  layer 

cv  consolidation  coefficient 

D  horizontal  dispersion  coefficient 

D50  median  floe  size 

Djj  rate  of  deformation  tensor 

Dm  settling  from  upper  layer 

E  total  wave  energy 

E '  wave  energy  density 

Eb  erosion  consolidated  bed 

Em  entrainment  rate 

ejjk  permutator:  e  =  1/-1  for  cyclical/anti-cyclical  /,y,  otherwise  e  =  0. 
Fj  wave-induced  traction  on  fluid  mud  (streaming) 

/  Coriolis  parameter 

fws  friction  coefficient  for  waves 

g  gravity 

H  wave  height 

h  thickness  of  upper  layer  («  water  depth) 
ho  initial  thickness  of  upper  layer 

Ko  coefficient  of  lateral  stress 

Kk  permeability  coefficient 

Kp  effective  stress  coefficient 

k  permeability 

k  wave  number 

kj  imaginary  part  wave  number  (wave  damping) 

k,  real  part  wave  number 

kx  coefficient  in  friction  coefficient 

kj  coefficient  in  erosion  formula 

k:  coefficient  in  formula  for  critical  shear  stress 

Me  erosion  parameter 

N  wave  action  density 

N  number  of  wave  cycles  during  liquefaction 

N{  number  of  wave  cycles  at  liquefaction 

n f  fractal  dimension 

PI  plasticity  index 

P,x  wave-induced  stresses  on  the  seabed 
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p  pressure 

p  isotropic  stress 

po  initial  isotropic  stress  (i.e.  prior  to  liquefaction) 

pw  amplitude  wave-induced  stresses  on  the  seabed 

Q  deviatoric  stress  integrated  over  fluid  mud  layer  thickness 

q  deviatoric  stress 

q  amplitude  deviatoric  stress 

R  complex  relative  wave  number 

Re  Reynolds  number 

Ri*  bulk  Richardson  number 

S  salinity 

S  switch  function:  S(.v)  =  .v  for  jc  >  0 
S'b  m  fluid  mud  induced  wave  dissipation 

Ty  stress  tensor 

t  time 

Uj  depth-mean  velocity  in  upper  layer 

u  velocity  of  upper  layer  in  x-direction 

He  Eulerian  velocity 

Ul  Glm  velocity  (Generalized  Lagrangean  Mean) 

Us  stokes  drift 

uorb  orbital  wave  motion 

u*  shear  velocity  of  main  flow 
w*/  shear  velocity  by  main  flow 
u*,m  shear  velocity  at  water-mud  interface 
u shear  velocity  by  waves 
um  velocity  of  fluid  mud  layer  in  x-direction 
v  flow  velocity  in  y-direction 

v,„  velocity  of  fluid  mud  layer  in  y-direction 
Wc  consolidation  velocity 

IVS  settling  velocity 

W  vertical  velocity 

w5  vertical  velocity  of  solids 

x  longitudinal  co-ordinate 

y  longitudinal  co-ordinate 

Zb  elevation  of  consolidated  bed  with  respect  to  reference  level 
Z,„  elevation  of  water-mud  interface  with  respect  to  reference  level 
Zs  elevation  of  water  surface  with  respect  to  reference  level 
I  Id  second  invariant  of  rate  of  deformation  tensor 

a  calibration  coefficient  fluid  mud  induced  dissipation 
8  thickness  fluid  mud  layer 

8mm  maximal  (user-defined)  thickness  fluid  mud  layer 
8)  initial  thickness  fluid  mud  layer 

A p  density  difference  over  water-mud  interface 

r r  turbulent  eddy  diffusivity 

y  shear  rate 

/)/  elevation  of  water  surface 

rj2  elevation  of  interface 
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6  wave  direction 

A  friction  coefficient 

Ag  complex  Stokes  number 

As  structure  parameter 

jU  mud  viscosity 

p '  mud  viscosity  parameter 

//x  residual  mud  viscosity  at  fully  deflocculated  state 
v  kinematic  viscosity 

vm  viscosity  mud  layer 

Vf  turbulent  eddy  viscosity 

p  fluid  density  of  upper  layer  (taking  into  account  the  effect  of  salinity) 

Pd,y  dry  bed  density  (mass  concentration)  of  fluid  mud  layer 
pm  bulk  density  of  fluid  mud  layer 

cr  wave  number  (frequency) 

cr  normal  and/or  principal  stress 

cc  vertical  effective  stress 

T  shear  stress 

Tb  bed  shear  stress 

Tb  Bingham  strength 

t’b  Bingham  strength  parameter 

rs.»  residual  Bingham  strength  at  fully  deflocculated  state 
Tc,e  critical  shear  stress  for  erosion 

r(/  critical  shear  stress  for  deposition 

Tz  shear  stress  at  level  Z  (i.e.  water  surface,  interface  and  base  of  mud  layer) 
v  Poison  ratio 

(/)  angel  of  repose  /  internal  friction 

<j)  volumetric  concentration 

(/)  phase  angle  water  elevation  and  flow  velocity  of  surface  waves 
</>o  angel  of  repose  /  internal  friction  at  fully  flocculated  state 
(j)'  phase  angle  between  surface  and  internal  waves 
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description 

wave  semi-orbital  excursion 

aggregation  coefficient  in  rate  of  flocculation  equation 

break-up  coefficient  in  rate  of  flocculation  equation 

sediment  concentration  by  mass  in  upper  layer 

cohesion 

cohesion  at  fully  flocculated  conditions 

wave  celerity  surface  wave 

wave  celerity  fluid  mud  -  water  interface 

group  velocity 

sediment  concentration  by  mass  in  fluid  mud  layer 

consolidation  coefficient 

horizontal  dispersion  coefficient 

median  floe  size 

rate  of  deformation  tensor 

deposition  rate  from  upper  layer 

total  wave  energy 

wave  energy  density 

erosion  rate  consolidated  bed 

erosion  rate  consolidated  bed 

entrainment  rate 

permutator:  e  =  1/-1  for  cyclical/anti-cyclical  i,  j,  k,  otherwise  c  =  0. 

wave-induced  traction  on  fluid  mud  (streaming) 

Coriolis  parameter 
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friction  coefficient  for  flow 

description 

friction  coefficient  for  flow 

friction  coefficient  for  waves 

gravity 

wave  height 

thickness  of  upper  layer  («  water  depth) 

mean  thickness  of  upper  layer 

coefficient  of  lateral  stress 

permeability  coefficient 

effective  stress  coefficient 

permeability 

wave  number 

imaginary  part  wave  number  (wave  damping) 

real  part  wave  number 

coefficient  in  friction  coefficient 

coefficient  in  erosion  formula 

coefficient  in  formula  for  critical  shear  stress 

erosion  parameter 

erosion  parameter 

wave  action  density 

number  of  wave  cycles  during  liquefaction 

number  of  wave  cycles  at  liquefaction 

fractal  dimension 

plasticity  index 
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wave-induced  stresses  on  the  seabed 

pressure 

description 

isotropic  stress 

initial  isotropic  stress  (i.e.  prior  to  liquefaction) 

amplitude  wave-induced  stresses  on  the  seabed 

critical  cyclical  stress  ratio 

dcviatoric  stress  integrated  over  fluid  mud  layer  thickness 

deviatoric  stress 

amplitude  deviatoric  stress 

complex  relative  wave  number 

Reynolds  number 

bulk  Richardson  number 

salinity 

switch  function:  S(x)  =  x  for  x  >  0 

fluid  mud  induced  wave  dissipation 

stress  tensor 

time 

depth-mean  velocity  in  upper  layer 

orbital  wave  motion  at  bed 

velocity  of  upper  layer  in  A-direction 

Eulerian  velocity 

GLM  velocity  (Generalized  Lagrangean  Mean)  in  .v-direction 

stokes  drift 

shear  velocity 
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shear  velocity  by  main  flow  at  water-mud  interface 

shear  velocity  at  bed-fluid  mud  interface 

description 

shear  velocity  by  waves  at  water-mud  interface 

velocity  of  fluid  mud  layer  in  .v-direction 

flow  velocity  in  v-direction 

velocity  of  fluid  mud  layer  in  y  -direction 

GLM  velocity  (Generalized  Lagrangean  Mean)  in  y-direction 

consolidation  rate  lower  layer 

settling  velocity 

vertical  velocity 

vertical  velocity  of  solids 

longitudinal  co-ordinate 

longitudinal  co-ordinate 

elevation  of  consolidated  bed  with  respect  to  reference  level 

elevation  of  water-mud  interface  with  respect  to  reference  level 

elevation  of  water  surface  with  respect  to  reference  level 

second  invariant  of  rate  of  deformation  tensor 

calibration  coefficient  fluid  mud  induced  dissipation 

thickness  fluid  mud  layer 

maximal  (user-defined)  thickness  fluid  mud  layer 

mean  thickness  fluid  mud  layer 

density  difference  over  water-mud  interface 

turbulent  eddy  diffusivity 

shear  rate 

elevation  of  water  surface 
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lower  layer 

used  in 

lower  layer 

lower  layer 

rheology 

rheology, 

lower  layer 

rheology, 

lower  layer 

SWAN, 

lower  layer 

rheology. 

SWAN 

computed  in 
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rheology 
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value 
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amplitude  elevation  of  water  surface 

elevation  of  interface 

description 

elevation  of  water  surface 

wave  direction 

friction  coefficient 

complex  Stokes  number 

structure  parameter 

mud  viscosity 

mud  viscosity  parameter 

residual  mud  viscosity  at  fully  deflocculatcd  state 

kinematic  viscosity 

viscosity  mud  layer 

turbulent  eddy  viscosity 

fluid  density  of  upper  layer  (taking  into  account  the  effect  of  salinity) 

dry  bed  density  (mass  concentration)  of  fluid  mud  layer  or  consolidated  bed 

bulk  density  of  fluid  mud  layer 

bulk  density  of  fluid  mud  layer 

bulk  density  of  fluid  mud  layer 

wave  number  (frequency) 

normal  and/or  principal  stress 

vertical  effective  stress 

shear  stress 

bed  shear  stress 

Bingham  strength 
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Bingham  strength  parameter 

description 

residual  Bingham  strength  at  fully  deflocculatcd  state 

critical  shear  stress  for  erosion 

shear  stress  at  level  Z  (i.e.  water  surface,  interface  and  base  of  mud  layer) 

Poison  ratio 

volumetric  concentration 

volumetric  concentration  prior  to  swelling 

phase  angle  between  surface  and  internal  waves 

angel  of  repose  /  internal  friction 

angel  of  repose  /  internal  friction 

phase  angle  water  elevation  and  flow  velocity  of  surface  waves 
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